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1. INTRODUCTION 

The major goals of the extension of NASA Contract 
NAS 9-12319 were concerned with building two general purpose 
numerical integration schemes into the NASA- JSC computer 
system. There were other minor goals associated with the 
contract and these will be discussed in Chapters 3 and 6. 

Since most of the work was concerned with the numerical 
integration schemes, the major portion of this report will be 
devoted to describing the state-of-the-art of numerical 
integration, the particular integrators built into the JSC 
computer system, and the use of the new integration packages. 
Those who are only interested in determining how to use the 
integrators may proceed immediately to the Appendices, which 
are self-contained. ■ ‘ 

The remaining portion of the report is as follows: 

Chapter 2 contains comments about numerical integration in 
general; Chapter 3 discusses the extrapolation numerical 
integration technique; Chapter 4 discusses the variable- 
order, variable-stepsize Adams numerical integration technique; 
Chapter 3 presents results concerned with numerical integration 
and optimization in the JSC PEACE parameter optimization program 
and Chapter 6 presents conclusions and recommendations. 
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2. GENUAL COMMENTS ON miKSRICAL INTEGRATORS 

In tills charier we shall briefly discuss the state-of-the-art 
of numerical integration, especially with respect to the 
development of general purpose numerical integration subroutines. 
For additional information, References 1 and 2 should be consulted. 
2,1 General Purpose Integrators 

As with any engineering system, the selection of a 
numerical integration subroutine usually involves a tradeoff 
between stability and performance. However, in the past few 
years a number of stable (reliable), relatively high-performance 
general purpose numerical integration subroutines have been 
developed. It the purpose of this report to introduce these 
routines to the Mission Planning and Analysis Division in as. 
simple a format as possible so that the routines will be used. 

As noted in Chapter 1 , the user, may proceed directly to the 
Appendices if he only interested in learning how to use the 
subroutines. 

The subroutines of this report are referred to as general 
purpose subroutines because the schemes are relatively flexible 
and apply to a large class of problems. . The schemes derive 
their flexibility mainly from the fact that, they are both 
variable order and variable stepsize. This fact allows the 
routines to adapt efficiently to many different physical situa- 
tions. In addition 5 the numerical integration routine of 
Chapter A (and Appendix B) has excellent diagnostic capabilities 
which indicate nn^erical .problems, e.g., excessive roundoff errors 
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One of the first questions which must be answered is: 

"When should I use a general purpose integrator?" The 
integrators described in Chapters- 3 , and 4 are very reliable, 
especially DVDQ. (Chap. 4-), and' they have been used with ease on 
problems from many different disciplines. If a new problem 
requiring numerical integration is to be solved, a general 
purpose integrator will not only produce reliable answers in 
a short time (with a minimum of expensive and tedious human 
effort) but usually also teach the user something about the 
physical problem, e.g., the stepsise usually contracts in 
areas of high acceleration or "action." Thus, -it is advisable 
to use a general purpose integrator when solving a physical 
problem for the first time . 

As to whether or not one should continue to use a general, 
purpose integrator, after fully understanding the process of 
the physical situation, is problem dependent. That is, since 
a general purpose integrator is in some sense conservative - 
(stable), one could probably invent rules- of- thumb to be used 
with a problem-dependent integrator to give a high-performance 
integrator for that particular problem. (This, of course, is 
the case with problems in celestial mechanics where many 
efficient, problem dependent integration schemes are employed.) 
In this process one must alv/ays keep in mind the tradeoff of 
•human effort and computer effort, e.g., it would not be 
worthwhile to spend many man-hours inventing rules- of- thumb 
for a particular problem if the problem Is to be solved only a 
few times, whereas the opposite would be true if thousands of 
simulations are anticipated. 
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Z.Z Main Existing Procedures 

The numerical integration schemes which are most often 

used in engineering analyses are the various fixed stepsize 

fourth-order Runge-Kutta and Adams predictor-corrector 

methods, } especially the Runge-Kutta methods. One reason 

for this is that the Runge-Kutta method is easy to understand, 

easy to program, and relatively reliable (once a workable 

stepsize has been selected) . Comparable properties of these 

two approaches (for fourth-order only) are listed below; 

Runge-Kutta: Based on Taylor series expansion; 

self-starting; requires four function evaluations 
per step; functions evaluated on a symmetrical, 
but unequal grid; difficult to approximate the 
local truncation error, which makes automatic 
variable stepsizing difficult. 

Adams Predictor-Corrector : Based on interpolation 

formulas; not self-starting, usually requires 
fourth-order Runge-Kutta scheme to develop the 
first few steps; after the starting phase, requires 
two* function evaluations per step; functions 
evaluated on an equal grid; relatively ease to 
approximate the local truncation error, v/hich makes 
automatic variable stepsizing relatively easy; 

.usually less stable than the Runge-Kutta method. 

From the comparisions above, one can see that the 

Runge-Kutta (or single step) and predictor-corrector (or 

multistep) methods are basically different, and thus, the 

performance of an integrator is problem dependent. Also, at 

first glance, the predictor-corrector method may appear to be 

faster because it only requires half as many function evaluations 

(after startup) as the Runge-Kutta method. However, this is not 


*This number may be higher in some schemes when more 
than one correction is allowed. 


necessarily the case because on some problems a smaller 
step-size is required for the predictor-corrector (because it 
is less stable) and/or the "overhead” (time associated with 
carrying out the requirements of the integrator) of a predictor 
corrector scheme is higher than for a corresponding Runge-Kutta 
scheme. ■ 

With regard to a variable stepsize, which is a necessity 
in many problems, one usually uses physical reasoning or a 
halving -and -doubling procedure to adjust the stepsize in a 
Runge-Kutta scheme. However, in the past few years a more 
attractive technique has been developed, and this will be 
described in Chapter Z.Jt. By comparing predictor and corrector 
values in a predictor-corrector scheme, it is relatively easy 
to adjust the stepsize. . 

2.3 Recent Developments 

In the opinion of this author, the three main developments 
in numerical integration in the past eight years have been the 
following (in chronological order): • 

(i) The development of an efficient rational function 

extrapolation numerical integration scheme;-^ 

(ii) The development of higher-order Runge-Kutta 

formulas with relatively efficient means for stepsize 
u 5 

control, ’ 

(iii) The development of high-order,, variable-order, 

variable-stepsize Adams methods in user oriented 

1 6 

subroutine form. ’ 

This report is based on the advances noted in (i) and (iii) . 



The higher-order Runge-Kutta schemes, sometimes called 
Fehlberg integration, are not included because Z4PAD already 
has a capability in this area. 

All of the general purpose integrators use a variable 

stepsize (although they can be used with a fixed stepsize) . 

The basis for the adjustment of the stepsize in (i) and (iii) 

will be discussed in Chapters 3 &nd 4* The means for efficient 

k 3 

adjustment of stepsize in Runge-Kutta schemes^’^ is as follows. 

4 

Runge-Kutta formulas can be developed for any order of accuracy, 
and for a given order the formula is not unique. For example, 
there exist numerous well-known fourth-order Runge-Kutta 
formulas. It can be shown that the least number of function 
evaluations required for a Runge-Kutta formula of order four is 
four, of order five is six, of order six is seven, of order 
seven is nine, and so on. Of course, formulas exist for the 
orders mentioned above which require more than the least number 
mentioned, and this is the basis for the efficient estimation 
of the stepsize. For example, suppose we wish to develop a 
■fifth-order Runge-Kutta scheme with a reliable stepsize 
adjustment procedure, and we select a fifth-order formula which 
requires the least number of function evaluations possible, 
i.e., six. A sixth-order formula requires at least seven function 
evaluations. However, suppose we construct a sixth-order formula 
with parameter values which match those of the fifth-order 
formula on the same six function evaluations required by the 
fifth-order formula, i.e., we embed the fifth-order formula in 
: the sixth-order formula. If this is done, then the remaining 



terms in the sixth-order formula give an excellent estimate 
of the truncation error for the fifth-order formula, which 
can then he used to adjust the stepsize. Fehlberg^ has shown 
that the resultant sixth-order formula requires eight function 
evaluations (whereas the minimum number for an arbitrary 
sixth-order formula is seven) . In summary then, for a fifth- 
order Kunge-Kutta • formula with the Fehlberg stepsize modification 
procedure, eight function evaluations are required (as opposed 
to six for a fixed stepsize fifth-order formula) . 

To conclude this section, we mention some recent comparative 

7-12 

studies of numerical integration schemes. The first 

statement to be made is that it is very difficult to draw 
definite conclusions from these studies. For example, the 
extrapolation scheme will perform n below par T? with a very 
poor choice of initial stepsize estimate, and thus, this 
scheme performs very well in Refs. 7-10, 12, but relatively 
poorly in Ref. 11 ► Also, some reports use the specified error 
tolerance as the accuracy indicator whereas . others compare the 
number of significant digits (compared to either an analytical 
solution or a finely-tuned numerical solution which is assumed 
to be exact to a large number of digits) . These two methods of 
comparison can cause differing conclusions. 

Another important point in comparing integrators is the 
.accuracy required. For low accuracies, e.g., 3 to 4 significant 
digits, the sophisticated integration packages will probably 
perform poorly. . However, in aerospace applications, high 
accuracy is usually desired and the situation reverses. 



Finally, a summary of major considerations in selecting 
or comparing integration schemes is listed below, assuming 
comparable accuracies for the. integrators . 

(1) CPU Time : This is probably the most important 
parameter to MPAD users since it indicates ho w fast the 
integrator is, i.e., how much computer time is involved. 

(2) Number of Function Evaluations : This parameter 

indicates how many times the right-hand sides of the 
differential equations are evaluated in a given run. 

In many comparisons this parameter is given the most 
emphasis (as opposed to CPU time) with the assumption 
that the smaller the number, the better the integrator. 
However, the scheme with the least number of function 
evaluations may not be the fastest (i.e., least CPU time) 
because of more overhead. The importance of the number 
of function evaluations parameter increases as the 
complexity of the right-hand sides of the differential 
equations increases. For example, an integrator v/ith 
libbl@ overhead, large number of function evaluations 
may be best when a spherical earth is assumed in a 
gravitational model, whereas an integrator with large 
overhead, small number of function evaluations is probably 
best if an extensive oblate model is employed. 

(3) Overhead : It is hard to put a number on this 

property of an integrator because it has to do v/ith the 
time required to carry-out the logical structure of the 
integrator, e.g., checks, rules- of- thumb , etc. built into 



the integrator. Usually, Runge-Kutta methods involve 
less overhead than predictor-corrector methods, and the 
DVDQ (Chap. 4) integration package involves more overhead 
than any other method reported in the comparisons in 
Refs. 7-12. However, on many problems it is the fastest 
integrator because of the small number of function 
evaluations required, which is due to the extensive 
logical structure built into the program. 

(4) Storage : In some cases the amount of computer 

storage required for an integrator is of importance. As 
integrators are developed which are applicable to almost 
all problems, this may be of some concern for small computer 
systems (e.g., DVDQ requires much more storage than other 
popular integrators; see Appendix B) - 

( 3 ) Scaling : The scaling of a physical problem can make 

or break some integrators. We have found that the routines 
of Chapters, 3 and 4 retain many of their desirable 
properties even on poorly scaled problems, which is not 
the case for most other integrators. 

( 6 ) Stiffness : A subject of considerable interest in 

numerical integration in recent years is "stiff" systems 
of differential equations. A stiff system is one which 
possesses at least two variables whose natural time scales 

■Z 

are significantly different (on the order of ICr), e.g,, a 
linear system with characteristic roots of - 1 , - 1000 . 

See Refs, 1 and 2 for further discussions of stiff systems. 



Reference 1 contains a subroutine for such systems, and 

F. T. Krogh of the Jet Propulsion Laboratory is currently 

1 3 

modifying DVLQ to include stiff systems. 



3- THE EXTRAPOLATION METHOD 
The goal of this chapter and Chapter if is not to develop 

4 f 

the underlying equations in the integration packages since 
these are readily available in textbook form, e.g., Refs. 1 
and 2. Instead we shall mainly discuss the basis for the 
methods and present simple examples to illustrate the main 
features. It is assumed that the reader is already familiar 
with the basic concepts of numerical integration, e.g., knowledge 
of the workings of the basic Runge— Kutta and predictor-corrector 
schemes . 

5.1 History and Basic Motivation of the Method 

The simplest formula for the numerical integration of the 
scalar differential equation 

X = f(t,x) , x(t Q )r:X o , , (3-D 

is. Euler's formula . •. 

x(t k+1 ) = x(t k ) + hf(t k ,x k ), (k=0j1>* • *) .(3-2) 

where h is the stepsise and x(t^) denotes the value of x obtained 
at the k— step in the process. This is a first-order formula, 
ana it can be shown 1 that the global error (i.e., the difference 
between the true and approximate solutions) can be represented 

by 

error at t k =x T (t k )-x h (t k )=:hE(t k )+0(h 2 ) , (3-3) 

T 

where x (tj ) denotes the true solution, x k (t k ) is the numerical 
solution for stepsize h, and E(t k ) is a portion of the error 
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which is a. function of the differential equation only if 
round-off error is neglected. In 1927, Richardson 1 ^ noted the 
following property: suppose (3.2) is used twice, first with 

a stepsize of h and second with a stepsize of h/2. Then, the 
©jrj’Op £op ©sell cs.lcu.Xs. tiou must 1© of Hi© following fopm* 

xT(t k ) " x h ( t k )=hE(t k ) + 0<h 2 ) (3-4) 

xT( V - V2 (t k ) = (h/2 ) E (V +0(h2) - (3.5) 

f 

Elimination of E(t^) from these two equations results in 

xT(t k )=2 V2 (t k ) “ x h (t k )+0(h2 )> (3-6) 

i.e., a second-order approximation has been obtained by combining 
two first-order approximations. In many cases such a procedure 
gives higher accuracy with fewer function evaluations. For 


example, consider the integration of x=-x from t =0 to t =1 - 

... . 0 £ 
vith x(0) = 1. Hie results using Euler's formula (i.e., Eq. (3.2)) 

and the "Richardson deferred-approach-to-the-limit" or the 


"extrapolation" approach (e.g., 


(3.6)) are shown in Table 3.1. 


h 


1 

1/2 

1/4 

1/8 

1/16 

1/32 

1/64 


x k (l) Euler's Formula 
($ of fn. eval.) 

0.0000000 ( 1 ) 
0.2500000 ( 2 ) 
0.3164063 (4) 
0.3435089 (8) 

0.3560741 (16) 
0.3620552 (32) 
0.3649865 (64) 


x h (1) Richardson 
(# of fn. eval.) 


0.5000000 (3) 
0.3828125 (7) 
0.3708116 (15) 
0.368 5393 (31) 
0.3680364 (63) 
0.3679177 (12?) 


TABLn 3.1 Comparison of Euler and Richardson's Deferred-Approach- 
to-the-Limit for x=~x. 


The true solution to the above problem is x(1 ) ^0 .367879 
(correct to six places). A simple computation shows that 
the extrapolated value corresponding to fifteen function 
evaluations is comparable in accuracy to the Euler value 
requiring sixty-four function evaluations. Thus, in this 
case quite a gain is achieved with the Richardson idea of 
. extrapolating values from results with a decreasing sequence 
of stepsizes. 

Although the above idea appears very attractive, it was 
somewhat dormant until W. B. Gragg studied the method 
extensively in his doctoral dissertation 1 ^ in 7 963 * The delay 
in the use of the idea was probably mainly due to the lack of 
a digital computer. Then, based upon the theoretical studies 
of Gragg, Bulirsch and Stoer^ developed an efficient numerical 
integration scheme involving this idea. The routine presented 
in Appendix A is a recently improved version of the Ref. 3 
subroutine. Since the scheme is relatively new, compared to 
Runge-Kutta and predictor-corrector methods, a good deal of 
current research in numerical analysis is concerned with the 
understanding and application of the extrapolation method. TIius, 
the scheme should be improved even more in the near future. 

3 . 2 . Motivation for Order Increasing Capability 

As the example of the previous section demonstrated, the 
order of approximation was improved by one when the extrapolation 
was performed on two integrations of the first-order Euler 
formula. Actually each additional integration of Euler’s 
method, with a monotonically decreasing stepsize, followed by 


extrapolation increases the order of approximation by one, as 
in the generation of Table 3A. The reason for this is because 
the error term for Euler’s method possesses an asymptotic 
expansion of the form 

x(t p ;h)ix T (t p )+E 1 (t p )h+El,(t p )h 2 +...+E m (t p )h m +0(h m+1 ), (3. 

where x(t ;h) is the result of the Euler integration at t with 
P P 

T f \ 

stepsize h, x (t ) is the true solution at t , and E.(t ), ... , 

P P 1 P 

E (t ) are portions of the error which depend only upon the 
m p 

differential equation (and not h.) « In Eq. (3.7) > the quantities 

- rn 

x (t ),E^(t ), . . . , E m (t ) are ? in general, .unknown. 

Consider the use of Euler’s formula on m+1 integrations 
with stepsizes h^> h^ > . . . >h .• Then, neglecting terms of 
order 0(h ra+ ^), Eq. (3.7) could be evaluated with each of the 
stepsizes (where the left-hand side of the equation is the value - 
of x (t^) determined by the integrations) , and the result is 
m+1 linear equations in the unknowns x (t^),E 1 (t ),..., E (t ) . 

• p I Tj lii U 

T th . - - . 

The resultant value for x (t ) is an m — -order approximation 

of the true solution. 

In practice one need not solve the system of linear equations 
mentioned above to obtain the higher-order approximation of 
x (t ). The sole reason for presenting the above example was 
to show how successive integrations increase the order by one 
if the Euler formula is the base formula. 

Of course one may use the extrapolation idea with any 
integration formula as the basic formula. Thus, the question 
arises as to the existence of a best base formula. The answer 
to such a question would involve to some degree finding a 



formula which increases the order of approximation by more than 
one with each integration; that is, the formula would possess 
an error terra with an asymptotic expansion of the form: 

X (t p ;h)=x T (t p )* |l r E^.(t p )h c ' k . . (c > 1 ) (3-8) 

Gragg 1 5 showed that if the modified mid-point rule is used as 
the basic formula, then r=1 , c =2 in Sq. : (3-8)> i.e., the error 
for the modified raid-point rule has an asymptotic expansion 

of the form: 

(3-9) 


x(t p ;h) = x T (t p ) + g- t p )h 2k . 

Gragg 1 s 1 ^ modified midpoint rule is as follows: 

h -H/n. , JT r , 7 - set of even increasing integers 

i i L 1 J 


X 


x ^ t o ) 


( 3 . 10 ) 


; v - x i = v h i f(t o’ x o ) • - ; 

Vs = V'V'Vi’Vi' ’ p=0, • • • ’ ni_1 

. x(t„+H;h.) = 0/h)x n +1 + 0/2)x n +OA)x i , 

U . -t- 1 ^ 

where H is the basic stepsize (e.g., the output interval). The 
integers in { n^ must he either all even or all odd for 
theoretical reasons; Gragg has shown that an even set has certain 
advantages over an odd set, and typical choices are { 2 , if, 6, 8, 

3 

12,16,21*,...} and £2,4,8, 1 6,32,64, •• •!• 

Finally, one other theoretical result due to Gragg is that 
the error between the true solution and the last extrapolated 
value, say x£(t ), is of the form 

K<V* T M S (3 '" ) 



where h > h*> ...>h , if the modified midpoint formula is 
o 1 m 

employed . 

3.5 Implementation Notes 

The heart of the extrapolation numerical integration 
scheme is a table of values which we shall call the "T- table", 
after Bulirsch and Stoer.^ To understand the operation of the 
method it is probably best to consider an example of its operation 
Example : As with our previous examples, assume that the base 
formula is Euler 1 s method. In the previous section we showed 
that an m^ order approximation could be obtained by solving 
a system of m+1 linear equations. Actually, the solution of the 
linear system can be avoided by employing a recursive extrapola- 
tion process due to Aitken.^ Let us first write the resultant 
table and recursive formula, and then use a simple example to 


motivate the result. ' 

Suppose four integrations with stepsizes V h 1 > h 2 > h 3 
are employed on the interval [t ^,.t^] with the resultant 
values x(t ;h Q ), x(t x(t ;h 2 ), x(t ? ;h^) . The following 

T-table is then constructed by the recursive formula 

4 - 4H ♦ ( 5.i2) 
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FIGURE 3 • 1 • T~ Table for Polynomial Extrapolation with 
Euler 1 s Method. 

In the table above only the first column is computed by 
numerical integration. The remaining columns are generated 
by the algebraic equation (3.12). Thus, with four integrations 
of the first-order Euler f s method, a fourth-order result, % T^, 
can be inferred. In fact, if one sets m^3 in Eq. (3-7) and 
evaluates the equation at h=h Q ,h^ , then the solution of 
the four equations for x (t ) will be precisely T^, which is 
obtained by a well-defined recursive process. Intuitively 
one would expect that the solution becomes more accurate as 
one either moves down a column or moves to the right along a 
diagonal. Gragg ^ has proved that under reasonable conditions, 
the T- table for polynomial extrapolation with the modified 
midpoint rule has the following properties: 

(i) The order n column converges to the solution 
faster than the order m column, with n>m. 



(ii) The principal diagonal, i.e., , . . . , T^ , . . . , 

converges to the solution faster than any column. 

Thus, one should take the "bottom element from the last column 
as the "best estimate of the solution, e.g., T° in the example 
of Figure 3 * 1 . 

To fix the interpretation of the T~ table, let us perform 
some simple computations. We shall show why T^ is a third- 
order approximation, and verify that T^ obtained by the recursive 
relation is the same as the value obtained by solving the 
corresponding system of linear equations. 

The minimum number of . integrations for a third-order 
approximation with Euler f s method is three. Thus, assuming 
h Q =h,h<j =h/2,h2=h/4 we have (at t ) : 

x(h) =x T +E 1 h+E 2 h 2 +0 (lib 

x(h/2)=x T +E 1 h/2+E 2 h 2 A+0(h 3 ) (3 - 1 3 > 

x(hA) =x T +E 1 h/if+E 2 h 2 /1 6+0 ( h 3 ) 

One can easily verify that Eq. (3. IE) implies •' • ".v: •. . 

T°=2x(h/2)-x(h) - . • ' \Y :: ■ 

t| = 2x(h/A)-x(h/2) , (3»14) 

which correspond to Eq. (3.6) developed earlier. An alternate 
interpretation of Eqs. (3.1^) is that they represent Eq. (3.13) 

p 

with E^ eliminated and 0(h ) terms neglected. In fact the 
elimination of E^ from Eqs. (3-13) implies 
x T =2x(h/2)-x(h)+h 2 E 2 /2+0(h 3 ) 

x T =2x(hA)-x(h/2)+h 2 E 2 /8+0(h 3 ). (3-13) 

Then, upon elimination of E^ from Eqs. (3-15 ) > ihe following 
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third-order approximation of x is obtained 

x T =[x(h)-6x(h/2)+8x(hA)]/3+0(h 3 ). (3-1 

This value, of course, corresponds to the value obtained 
by repeated use of Eq. (3*12). 

Before v/e consider variable stepsise and order aspects of 
the method, note that the name ’’polynomial extrapolation” is 
due to the fact that the value of x(t ;0) is obtained by using 
a polynomial fit of the data x( t jh ) , x(t jh* ) , ... to 

jr ^ ir 1 

extrapolate the value at h=0. It is not necessary to use a 
polynomial fit of the data, and in fact, the Bulirsch-Stoer 
extrapolation scheme uses a rational function fit of the data. 
In most simulations to date, the rational function procedure 
has given better results. In the case of rational function 
extrapolation v/ith the midpoint rule, the form of the T- table 
remains the same, however the recursive formula (3.12) is 
replaced by^ 

T^x(t ;\) (3. 


mi rpi + 1 , 

Wi 


mi**" 1 n>l* 


( h i V 1- 

\ h i+k / 


rpi+l 

x k-1 


mi+1 

x k-2 


-1 


. (k> 1 ) 


These are the formulas in the subroutine of Appendix A. 



5.4 Variable -Stepsize , Variable-Order Considerations 


The main reason for using variable stepsize and order 

methods is to develop .the most efficient scheme possible for a 

specified accuracy. From the previous sections it is apparent 

that more than one stepsize is utilized in the basic operation 

of the extrapolation algorithm, and the only questions with 

respect to stepsize have to do with the initial choice and the 

final choice, e.g., h and h if [h /2 ,h /4 ? . - . >h =h /2m] is 

o m o o mo 

the sequence. 

The comparative studies in Kefs. 7“]] indicate that a 

poor choice for h^ can cause inefficiency in the scheme, 

e.g., if h Q is too large the extrapolated values may be 

contaminated by round-off and if h is too small it is not 

o 

talking full advantage of the extrapolation process. Both the 
scheme of Appendix A and the original version (pp. 96-99 ? 

Ref. 1) have means for adjusting h Q by checking how many values 
are required for the first column of the T- table (i.e., if too 
few are required, the stepsize is too small; if too many are 
required, the stepsize is too large) . 

The determination of order and the final stepsize, h , 
are somewhat coupled. The algorithm of Appendix A can generate 
six columns, which implies a possible twelth-order scheme since 
the orders of the columns go up by twos with the midpoint formula 
(see Eqs. (3*9), (3*10)). The termination of integrations 
(or the determination of h ) varies from basic step to basic 
step. Integration is terminated when two successive approxima- 
tions T^ - ^ and p^“* v+ 1 differ by less than the specified error 



tolerance, and the value is taken as the best estimate. 

This feature is one of the strongest aspects of the scheme 
since the cutoff is determined by comparing successive better 
approximations to the solution. One may consult the table 
on page 90 of Ref. 1 to see a striking example of this feature 



4. VARIABLE ORDER, VARIABLE 
STEPSIZE INTEGRATORS 


In this chapter we shall discuss the basic ideas behind 

variable order, variable stepsize Adams . predictor-corrector 

methods. Although there exist relatively efficient variable 

stepsize Runge-Kutta integrators, e.g., Ref. 4> there does not 

appear to be an efficient way of changing the order in a 

Runge-Kutta scheme. The two main variable order, variable 

& 1 

stepsize Adams methods are those due to Krogh and Gear . 
Krogh f s program is discussed and listed in Appendix B, while 
Gear's program is listed on pp. 158-166 of Ref. 1. 

4.1 Storage Methods 

Before we discuss means of storing information available 
in a predictor-corrector method, let us write the first few 
Adams formulas for future reference. -With 

* = f(t lX ), Xp =x(t p ),fp=(t p) x p ), h=t p+1 — t p (4.1) 

we have the following. 


Adams-Bashforth (Predictors) 

Error 

VrV h f P 

(h 2 / 2 )x 

x p+ 1 -x p = d 1 / 2 5 f p- fp_ •] ) 

(5h 3 / 12 )x 

Vi- x P =(h/12)(23 V 16 Vi +5f P-2 ) 

(9hV 24 )x (4) 

VrV (h/ 24 } ( 55 V 59f P-i +37f P-2- 9f P-3 ) 

( 251 h 5 / ?20 )x 
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Adams -Moult on (Correctors) 


Error 


Vi~v h Vi 

-(h 2 / 2 )x 

VrV (h/2)( Vi + V 

-(h 5 / 12 )x 


-(hV 2Zf )x (if) 

Vr V (h/2Zf) (9 V 1 + 1 9 V 5 Vi + VP 

• i 

• • . 

-(19h 5 / ?20 )x (5) 

• 

• ♦ 

ft. 

The above formulas are used with a k— (or k> 

4 

-1 ) order predictor 

t h Ti 

and k — order corrector, where a formula is k— order if its 

error term is of the form h k+1 

Classically, the 

Adams . predictor and corrector formulas v/ere 

l r 

written as 

x - + h ? ^ a .f 

p+1 p 3=0 3 p-3 

(Predictor) ( 4 - 2 ) 

k * 


+ h a. f . , 

P +1 P .3 p + 1 ~3 

(Corrector) ( 4 - 3 ) 


and stored in the computer in the form x ,f , f ...... -f , . 

P P p-1 3 p-k 

Actually one loses a great deal of free information by 'choosing 
this method for storing hack information. Let us first state 
two other means for storing hack information, and then discuss 
the advantages of using these methods. 

Backward Difference Table 

Store: x p , V 1-1 f p (4.4) 

(Used in DVDQf the notation will be explained below.) 

Scaled Derivatives 

Store: z =(x ,hx , . . . , Ow , ) X (< 1 .)]T 

P P P / ql / p J 

(Used by GeaP and popularized- by Hordsieck.V. 


(if. 5) 
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It can be shovm that well-defined linear transformations connect 
all three means of representation."* Vie shall see below, by means 
of an example, how one would construct such a linear transformation. 

The Adams predictor-corrector formulas written in 
backward difference form are 

x^|=x + h2L q. V J f (fe^ order) (4»6) 

P*1 P j = o J p 



=V 


k-1 * 

kJS q, 

3=0 J 



Vi 5 


(k+1 -order) (4*7) 


where the different order predictor 
employed since these are precisely 
The notation is as follows: 

v °V f P ■ v J ' +1 V v 


and corrector formulas are 
the formulas used in DVDQ. 


4 p _ 1 (j2 0) • (if-. 8) 


v P f p =i1 Vi jx S ] 


■jtl 


P+1 


= V^f 


p P+1 


- V 3 fJjiO). (4-9) 


Let us consider an example which clarifies the notation and 
demonstrates how one would construct a linear transformation 
between the various means of representation, e.g., between 
Eqs . (4*2), (4-3), (4.6), and (4-7). 

py ample : Consider the second-order Adams-Bashforth formula, i.e., 

( 4 . 10 ) 


.(P) 


x'*' - x^ = (h/2)(3f ■ ,). 

p+1 p p p 1 

With respect to Eq. ( 4 . 6 ), k=2, so in backward difference form 


Eq. (4*10) can be vnritten as 

1 


■(P) 

p+1 


= x + h 


3=0 


q* v J f. 


(4.n) 


or 



(4.12) 




(P) 


1 " *n + 11 [q O V °V *1 


P+1 ”p ” K "O v “P J -1 r ”p- 

We shall now employ the notation definitions of Eq. (4-8) to 
write Eq. (4.12) in terms of function evaluations only, and then 
determine q Q and so that Eq. (4-11) matches Eq. (4-10). By 
definition: ' 


V°f =f 
v P P 


vh p = v°f p - 


:p-i p p-i 
Upon substitution into Eq. (if .12), we have 


or 


: (P) 

p+1 


■CP) 


= v h[q oV q i ( Wi )] 


x p ;t= v hC( v«i )f P ^ 1 Vt 1; 

Comparison of Eqs. ( 4 . 10) and (4*13) implies 
V = 3/2 , q 1 = 1/2 , 


(if. 13) 


(if-lif) 

which is the desired linear transformation between the coefficients 
of Eqs. (4*2) and (4*6) for the second-order case. ‘ • 

Now that we have seen that the backward difference 
representation of the predictor-corrector formulas is no more 
difficult than the standard means of representation, the question 
arises as to what we gain out of such a representation. The 
answer is given by noting the following equations: 


Vf= f =x 
P P P 


'''rt-Vi*«VrV M i ■ M(i p-r i p )/< - h>1 


P P P 
V 2 f p ^ (h 2 / 2 )W 


hi 




(if. 13 ) 



Equation ( 4 . 15) has great significance since the major portion 
of the local truncation error can he represented by (see Gear, 

p. Ill) 

C k+ 1 h k+ 1 x£ k+1) + 0(h k+2 ), (4.16) 

i.e., the higher-order differences indicate the local truncation 
errors for various order integrators since they are proportional 
to the higher-order derivatives of x at t^. Thus, by employing 
this means of representation, or the scaled derivative representa- 
tion, one has the means for not only adjusting the stepsize but 
also for choosing the most efficient order integration formulas. 
4-2 Variable Order, Variable Stepsize Procedures 

In the previous section we saw that by storing back informa- 
tion in either backward difference or scaled derivative form, 

--we have the means for estimating -the local- truncation error at 
various orders with no extra computation. In this section we 
shall present a method due to Gear for adjusting .the stepsize 
and order automatically. • - - • 

Suppose that the integration has proceeded long enough to 
generate a sufficiently long "tail 15 of scaled derivatives or 
backward differences. Later we shall discuss problems associated 

with starting the algorithms. Suppose we are currently at t with 

P 

k— order integration formulas, stepsize h, and scalar truncation 
error parameter E. Use the k— order formulas and stepsize h 
to compute where t^ + 1 =t^+h. Erom our "tail" of backward 

differences or scaled derivatives we can easily compute the 
truncation error associated with the k^ order formula, and make 
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the following check: 

r .k+l„(k+1) ? F 

C k + 1 h Vi - 


( 4 - 17 ) 


If Inequality ( 4 * 1 ?) is satisfied, we accept the computed value 
for x(t -j); if the inequality is not satisfied, we must determine 
a smaller stepsize h to compute a value x(t^ + ^) which satisfies 
( 4 . 17 ) . In either case we proceed to the tests below. 

We now wish to either possibly increase the stepsize if 
(4*17) was satisfied or decrease the stepsize if (4. 17) was 
not satisfied. In addition we wish to check to see if the 
order should be increased or decreased. Define: 

New stepsize = zh, z scalar unknown. ( 4 * 18 ) 

If all quantities were known exactly and the k+1 -derivative 
of x^ + ^ remained constant,, then the optimal stepsize would be 
defined by (4.17) with the equality, i.e.. 


r /^k+1 (k+1) P 

C k+ i(zh) x p+1 = E. 

Then, solving for z we obtain 


= {Mc k+1 h k+ U^]j 


( 4 . 19 ) 


(4.20) 


(k+1 ) 

However, the quantities are not knovm exactly and x v 7 does 
not, in general, remain constant so a safety factor is included 
in Eq. ( 4 * 20 ), e.g., 1 


^ = iM E /fw k+lx Si 1)] 3 


(4.21) 


where the k subscript is attached to z to imply that this is 
the stepsize multiplier associated with the k— order integration 
formula. Since we also have backward differences which imply 
and x^j^ , we can determine the optimal stepsizes at the 



neighboring orders, i.e., at the k-1 and k+1 orders. Operating 
in exactly the same way as we did in forming Eq. (4.21), we 
obtain (with appropriate safety factors) 

JL 

2 k-1 = T3 { E/ ! c k h ^!j ]] k (Order k-1) (4-22) 

1 

z k+1 = l!7i:l E/fC k+2 hk+2x p+1 2):i ^ k+2 • (Order k+1)(/+. 23) 
The safety factors 1.2, 1 .3 } . and 1 .4 are chosen so that we have 
the following order of priority: no order change, lower the 

order, and increase the order. The order which produces the 
largest value of z in Eqs. (4*21), (4*22), (4*23) is selected 
as the integration formula order for the next integration. 

With the procedure above, an arbitrary value for h may 
result, and thus to employ the current "tail" of information, 
interpolation must be used to shift the information to the 
new stepsize. As shown by Gear"*, this may be done by means of 
a matrix multiplication. Also, because, of the overhead associated 
with implementing the tests and interpolation, the above 
procedure is not performed on each step in Gear's program; see 
page 137 of Ref. 1 for a discussion of hov/ the procedure is 
implemented. 

The same ideas form the basis, for changing order and stepsize 
in DVDQ. However, DVDQ only allows halving and doubling of 
the stepsize, which Krogh claims is more efficient when the 
order is to be varied. Also, DVDQ surveys the local truncation 
error at four orders to build more stability into the order 
changing procedure. 
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4*3 Self-Starting Adams Methods 

To conclude this chapter, we shall present an example 
to demonstrate how predictor-corrector methods may he used 
without employing alternate integration schemes to get them 
started, e.g., many existing fourth-order predictor-corrector 
subroutines employ a fourth-order Runge-Kutta scheme to 
generate the required starting values. Before we consider 
the example, let us consider an alternate method for estima- 
ting the major portion of the local truncation error in a 
predictor-corrector scheme. Suppose the predictor and 

corrector are-of the same order. Then, denoting the true 

T 

solution of the differential equation at t by x , it can be 

jr ir 

shown that. (Ref. 1, p. Ill): 

x< P ) = + C,h k+1 x^ +1 ' + 0(h k+2 ^ • ( 4 * 24 ) 

p p k p ' * * ' 


x (C) = x T + C, h k+1 x (k+1) + 0(h k+2) 
P P k p 


T 


(4-25) 


Upon elimination of x between these two equations, we obtain 


. h 1 ‘*V< k *»(c k -c lt ) . 0 (h>'* 2 ) , 


(4.26 ) 


( k+ 1 ) 

where everything is known to order k+1 except x^ ' . Thus, 
( k+ 1 ) 

solving for x^ : 


x^ k+1) « [x< P) - x£ C) ]/[h k+ 1 (C k -C k )]. 


(4-2?) 


( k+ 1 ) 

Then, given an estimate of , we can estimate the difference 

ir 

between the true and corrected solutions from Eq. (4*25) > i.e., 


jx (C) ~x T | S I C, h k+1 x (kl1) 

' K P 


(4-28) 
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where x^* + ^ is estimated from Eq. (if. 27). Upper and lower 
P 

limits, and E^, respective^, are usually specified to 
indicate when the stepsise should be halved or doubled, e.g., 


} X (C >-X T | 


> E^ Halve stepsize 


< E^ ^ Double stepsize 


(4.29) 


(4.30) 


Let us now use this method in an example which demonstrates a 
self-starting Adams' procedure. 

Example : Given x^f(t,x), t Q ,x(t o )=x o , and error constants 

E-^, E^, and E.' Also, an initial stepsize h Q must be specified 
and this value should be smaller than one would expect to result 
with a' fourth-order integration scheme. For sake of illustration, 
assume that only one correction is employed after each predictor 
step. Consider the first-order Adams predictor and corrector 


formulas. Then, 
„(p) .. „ ^ 


= x +h f 
0 00 


(4-3D 


= VV C V X 1 ] - 


- ( 4 . 32 ) 


We now check to see if the stepsize should be changed by using 

#* o 

Eq. ( 4 . 27 ) to estimate , and then checking 


E < 
n — 


■ x i 1 - \ 


\< |'i-(x 1 (P) - xj C) )U E u . 


(4-33) 


After maiming this test, we either halve, double, or keep the 
same stepsize. Let h be the resultant stepsize. We then form the 



(4-34) 


first components of the backward difference table 

We now have enough information to use the second-order 
Adams^ formulas. Again using a single correction v/e form: 

X CP) = x 1 + (h/2)(3f r f 0 ) (4.35) 

x 2.^ = x 1 +(h/2)[f 2 (t 2 ,x| P ^) + f^]. (4.36) 

Of course if v/e are storing backward differences, v/e would use 
the backward difference version of the predictor-corrector 
equations in the actual calculation (i.e., Eqs.- ( 4 . 6 ), (4*7)) > 
Given x^ , we again employ Eq. ( 4 . 27 ) to estimate x^, and 
check to see if the stepsize should be changed by 

— 1 4 c) - x 2 i - v ; (4 .37) 

After this test, we either halve, double, or keep the same 
stepsize, and then augment the difference table to form 

[x| c ' , V°f 2 . V r f z , v 2 f 2 3- • (4-38) 

One then proceeds in the same manner to the third-order formula, 
and eventually the difference table will possess sufficient 
information to allow a switch over to the automatic order and 
stepsize changing procedure described in the previous section* 

Of course, if one desires a fixed order predictor-corrector 
method, e.g., fourth-order, the above starting procedures would 
be used up to fourth-order and then sufficient information will 
exist for the use of the fourth-order formulas. 

In Appendix B a listing of DVDO is presented along with a 


sample program. 


5. RESULTS WITH NASA- JSC COMPUTER PROGRAMS 

In this section further results involving the NASA- JSC 
PEACE parameter optimization program will be discussed. Since 
most of the work in this area was done in the first portion of 
the contract (as opposed to the extension) , and Is thoroughly- 
documented in Ref. 18, we shall only discuss recent results. 

5.1 Parameter Optimization 

The main modification to the PEACE program involved the 
addition of the Pletcher method with a one-dimensional - search. - 
The method was simulated on the stage -and -half configuration 
used for the simulations in Ref. 18. It was found that the 
Fletcher method with a one-dimensional search did not perform 
as well, as either the DFP or Broyden algorithms (see Table 5-1) 
It appears that the terminal convergence properties are poor 
because the H-matrix probably becomes H contaminated" by the 
switching of the formulas for the H-matrix. Note that, from 
Table 5-1? the rate of convergence of the Fletcher method 
slowed considerably after the first time that the formulas 
were switched. That is, the test within Fletcher's method 
implied that the DFP formula should be used until the sixth 
iterate when the Broyden formula was used. Thereafter the 
convergence was very slow. Recommendations for use of the 
schemes will be given in Chapter 6. 
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Iteration 

Number 

COST 

DFP 

(All with 1- 
Broyden 

<D Search) 

Fletcher 

0 

1611. 

1611. 

1611 . 

1 

271.1 

271 .1 

271.1 

2 

' 13*73 

13.73 

13.73 

3 

2.792 

2.792 

2.792 


1.311 

1 .311 

1.311 

3 

.8897 

.8897 

:8897 

6 

.2636 

.2636 

.2636- 

- 7 

. 1 621 

.1621 

.2462 

8 

.1563 

• .1563 

.1603 

9 

• 1554 

.1554 

.1580- 

10 

.1461 

.1461 

' .1580 

: 12 

.1183 

.1183 

.1372 

■•■14 

.1046 

-1043 

.1548 

16 

.0972 

• 0967 

*1472 

18 

.0923 

‘ .0921 

,1470 

20 

.0888 

.0886 

■ .1413 


*(In Fletcher column, * indicates Broyden formula used; 
otherwise DFP formula used. The Broyden formula was 
also used on the 13 th and 19 th iterates.) 


TABLE 3-1 Comparison of DFP, Broyden, and Fletcher (all with 
1-D search) on Stage-And-Half Configuration. 



5.2 Numerical Integration 


The extrapolation numerical integration scheme of 
Appendix A has been built into the stage -and- half configuration 
version of the PEACE parameter optimization program. Since the 
program had its fqrmer integration scheme (fourth-order Runge- 
Kutta) built into the program (as opposed to being subroutined) , 
the extrapolation scheme will probably not perform as efficiently 
as possible. However, the resultant integrations should give 
some indication of its efficiency in a problem with a considerable 
number of discontinuities, and it should give some indication as 
to the accuracy of the Runge-Kutta integrations previously obtained. 
The computer program will be tested by NAS A- JSC personnel. 



6. SUMMARY AI03 EECOMMETOATIONS 


6 . 1 Summary 

Two general purpose numerical integration computer programs 
have been delivered (and checked out) to NASA- J SC personnel. 
These are the Bulirsch-Stoer extrapolation scheme (1972 version) 
and Krogh f s variable -order, variable-stepsize Adams method. 
User's guides and listings of the programs are presented in the 
Appendices. 

The PEACE parameter optimisation program was modified to 
include Fletcher's method with a one-dimensional search and the 
extrapolation integrator. In addition, W. F.' Powers presented 
seven lectures on optimal control and numerical integration to 
Mission Planning and Analysis personnel in August 1973* 

6.2 Recommendations 

1.) With respect to parameter optimization, the Broyden 
and DFP algorithms are recommend_ed if good terminal 
convergence is desired. In this respect, Broyden' s 
method has always performed better than DFP, but 
not appreciably better. Fletcher's method (without 
a 1-D search) appears to work well in the early 
stages, and especially on problems where the H-matrix 
in the DFP method is having trouble. However, 
because of the H-matrix switching it appears to have 
trouble obtaining rapid terminal convergence. 
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The Broyden method appears to do naturally and 
continuously what the Fletcher method does roughly, 
i.e., take tendency to singularity ox the H-matrix 
into account. From a theoretical point of view, 
the Broyden method appears to he the most attractive. 
With respect to numerical integration, the DVDQ 
scheme has the hest set of diagnostics and it is 
strongly recommended for use on new problems , The 
resultant output should indicate to the user a 
natural partitioning of the problem since the 
changing stepsize and order- usually indicate a change 
in physical characteristics of the solution, also. 

One may then use this information to define fixed- 
stepsise, fixed-order formulas in the various phases 
of the physical problem if a large number of production 
runs are anticipated. That is, there is no use paying 
the overhead of DVDQ if all of the runs will be 
similar with respect to physical characteristics. 

The extrapolation numerical integration subroutine 
has less diagnostic capability than DVDQ, but appears 
to be faster on runs where the right-hand side of the 
differential equation is not unduly complicated. This 
scheme is relatively new and should be continuously 
improved by research in the next decade. It may be 
more optimal for MPAD problems than DVDQ because it 
is a one-step method which should not be as adversely 
affected by multi-stage problems as DVDQ:. 



It is not claimed that the numerical integration 
routines listed in the Appendices will speed up 
existing programs, which usually employ a fixed 
stepsize that changes from one v/ell-de fined physical 
phase of the problem to another. However, they 
should speed up the process of choosing optimal 
stepsizes and determining accurate reference cases 
which will aid in the development of production 
programs. In this respect, the possibility of 
including both routines in SVDS should be considered. 
Hough rule s-of- thumb for choosing a numerical 
integration scheme are the following: 

(i) if low accuracy is all that is required 
(e.g., three significant digits) and/or there 
exists a large number of physical phases (which 
require starting and stopping the integrator), 
then a low-order Runge-Kutta scheme (e.g.,. order 
one to four) with sufficiently small fixed 
stepsize is probably optimal; 

(ii) if high accuracy is required (e.g., five or 

more significant digits), then either DVDQ or 

extrapolation should be employed. The relative 

efficiency of these schemes increases with the 

order of accuracy required. To date, comparative 
7-12 

studies have shown the two schemes to have 
the following properties: DVDQ has higher over- 

■ head, but smaller number of function evaluations 



than extrapolation, and for many problems, 
roughly the same computer time. In Ref. 9 a. 
rule for choosing between the two is proposed: 
if the right-hand side of the differential 
equation is relatively simple, extrapolation 
will probably give the least CPU time, whereas 
DVDQ will give the least if the right-hand side 
is lengthy (e.g., if gravitational anomalies 
are taken into account) . 
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APPENDIX A 
DIFSTS 

USER'S GUIDE AND LISTING 

SUBROUTINE DIFSTS is a double precision rational function 
extrapolation numerical integration scheme. It is an improved 
version of the ALGOL subroutine reported in Ref. 3* The 
original subroutine is also documented in FORTRAN in Ref. 1 
(pp 96-99) • The version presented in this appendix was supplied 
by its developer, R. Bulirsch. Some user's may prefer to use 
the original version (presented in Ref. 1) because -it contains 
more comment cards and error checks. However, to date, we have 
not encountered any difficulties with the faster, more stream- 
lined version presented in this appendix. 

DESCRIPTION 

DIFSTS is called by: 

CALL DIFSTS (N, TF, EPS, H, X, T) • 

where : 

N=order of system (number of differential equations) 

TF=a user supplied subroutine which calculates the 
derivatives, and has the form 

SUBROUTINE YF(X, Y, DT) 

X = independent variable 

v s dependent variable vector (must be dimensioned 
in TF) 

DY = derivatives (right-hand side of the differential 
. equation vector; must be dimensioned in YF) 

(Note: YF must be declared in an EXTERNAL statement 

in the program which calls DIFSYS . ) 

A1 



EPS r= stepsize error control (DIFSYS will reset EPS to 
10~ if the user supplies a- smaller number) 

H t= maximum integration interval 

X = independent variable ' • 

Y- = dependent variable vector (must be dimensioned N in 
calling program) 

The' quantities N, EPS, H, X, and Y(N) must be supplied before 
DIFSYS is called. 

OPERATION . * 

This subroutine does only one integration step per call. 
Hence, if the interval of integration is [x ,x^], where x^ may 
be defined implicitly, the user must test for x^ (or the implicit 
condition) and adjust H at the end so that is satisfied 
exactly. 

EXAMPLE PROGRAM 

Consider the integration of: 

•• *2 . 

X e -X + tX 

t ,t~,x(t ),x(t ) specified, 
o f o o 

Define y^ = x, Then, the following program will execute 

the integration with DIFSYS. 



IMPLICIT RSAL*8(A-H,0-Z) 

EXTERNAL YF . 

DIMENSION Y (2) 

READ (5j 100) EPS, H, TO, TF , (Y(I) ,1=1 ,2) 

N=2 

T=TO 

CALL’ DIFSYS (N, YF,EPS,H,T, Y) 

WRITE(6, 101 )T, (Y(I) ,1=1 ,2) 
IF(T.GE.TF)S10P 
IF((T+H) ,GT.TF)H=TF-T 
GO TO 1 

100 FORMAT(6D13.6) 

101 FORMAT ( ' TIME= ' , D20 . 8 , 1 Y(I) = ' , 2D20 .8) 
END 

SUBROUTINE YF(T, Y,DY) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION Y(2) ,DY(2) 

DY(1 )=Y(2) 

DY(2)=-Y(2)**2+Y(1 )*T 
RETURN 


END 
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SUfH OUT I ML: PI F S Y S ( N , Y F , EPS ,H »X ,Y ) 

REAL *8 Y(N),YA(1U),YL(1C) , Yi V H 10) »D Y ( 9 ) t UZ(10) ,DT (10,7) 
1 0(7),S(I0 ),X,XN ? H t G,P»B 1 »U,V,CtTA,W 

REAL *4 EP ( 4 ) / 0 - 4£ - 1 , C - 166-2 , G*64E-4 , 0 . 2 56 £-5/ 

LOGICAL-*! KONV f BG,KL,GR 
JT I =0 
F Y= 1 - 

ET A= AGS ( EPS ) 

IF (E TA.LT* 1 . £-1 1 ) ETA=1.E-11 

DO 100 1 = 1, N 
100 YA(I) = Y{ I ) 

CALL YF ( X, Y,DZ) 

10 XN= X +H 

B0=. FALSE. 

DC 110 1=1, N 
110 S( I ) = 0.00 
. M=1 
JR = 2 
JS= 3 

DO 260 J=l , 10 
IF(.UOTcbO) GOTO 200 
D( 2 )=1 „ 77777777 777 7776DC 
D( 4) =7. 11 11 11 11 111 111 IDO 
D ( 6 ) = 2 * 8444444444 4444 4D 1 
GO TC 201 

200 DC2)=2. 2500 
0(4)=9.D0 

D ( 6 ) =3 * 6D1 

201 IF C J.LE.,7) GOTO 202 

' L=7 • , ' 

DC 7 ) =6 . 4D1 
GO TO 203 

202 L=J ' ' 

D(L) =M*M . . • 

203 KONV =L.G T.3 
M=M + M 
G=H/P 

8 = G+ G 

DO 210 1 = 1, N 
YL( I ) =Y A ( I) 

2ic ym< i )=ya( n+G*nzm 

M = N-1 

00 220 K = 1 , M 

CALL YF{ X+K*G t YM,DY) 

DO 220 1=1, N 
U=YL<I)+8*DY(I> 

YL ( I ) = YAH I ) 

YM( I )= U 
U=OA BS ( U ) 

IFtU-GT.S ill) S ( I ) = U 
220 CONTINUE 

CALL YF (XN,YM, DY) ■ 

K L = L • LT *2 
GR = L * GT« 5 
FS - 0 • 
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DO 2 33 I = 1 » N 
V -DT (M) 

C= ( YM( I) *YL ( I MG*DY ( I ) MO .51)0 
DT ( 1*1) =C 
TA = C 

IF{ KL) GC TC 233 
DC 231 K - 2 » L . - . . 

B 1=0 { K MV 
D=B1-C 

u' = C — V . . 

U=V 

IFCB.EQ.O.DO) GO TO 230 
8=W/B 
u=c*n 
C = S1*B 

23C V = DT ( I T K ) 

DT t I t K) =U 

231 TA=U+TA 

IF ( •NOT .KONV) GO TO 232 

IF(OABS(Y(r )-TA) „GT.S (IMETA) KCNV=. FALSE. 

232 IF( GR .OR ,S( I ) oEQ.O.DO) GO TO 233 
F V= 0 A8S (W)/S(I J 

IF C F 5 . LT . FV ) FS=FV • 

233 Y ( I )=TA * 

TFCFS .EQ.O.DO) GO TO 250 - - 

F A = F Y 
K = L-1 

FY= (EP(K)/FS M* (1 ./(UK)) 

IF ( L «. c 0 » 2 I GO TO 240 

I F( FY.LT *0 .7-FA ) GC TC 250 V 

24C- IF(FY.GT.0.7) GOTO 250 ‘ 

H=H*FY 

JTI=JTI+1 ' • * 

,IF( JTI.GT.5) GO TO 30 
GO TO 10 

25 C IF (KONV) GO TO 20 

DC 3 ) =4. DO 
D< 5 M1.6D1 
80= .NOT *8 C 

m=jr • ' 

JR = JS 

26 0 J S=M+M 

H = H* C • 5D 0 
GO TO 10 

20 X = XN . • . 

H=H*FY 

RETURN •' 

30 H = 0 • DO 

DO 300 1 = 1 1 N 
300 Y(l)=YA(l) 

RETURN 
EN D 



APPENDIX B 


DVDQ- 

USEE'S GUIDE AND LISTING 

SUBROUTINE DVDQ is a double precision variable order, variable 
stepsize Adams predictor-corrector numerical integration 

.developed by I. T. Krogh of the Jet Propulsion Laboratory. 
As’ one may see by inspecting the listing, the program is 
extremely well-documented and over half of the listing is devoted 
to comment cards. Before we discuss the implementation of the 
deck, some .notable features of the program are listed below: 

(i) Maximum integration formula order =16. 

(ii) Halves and doubles for variable stepsize. 

(iii) Only MAIN and DVDQ are necessary, i.e M one 
need not employ a separate subroutine for the right- 
hand sides of the differential .equations. 

(iv) Order of the predictor = order of the corrector -1 . 

(Recall the general property that 

x p=* ( p )+ 0(h k+1 ) 

Xp = + 0(h k+m+1 ) + 0(h r+1 ) 

where k=order of predictor, r = order of the corrector, 
m = number of applications of the corrector. This 
implies that with k = r or r-1 the same order of 
accuracy is obtained for a single corrector application.) 

(v) Has "GSTOP" feature, i.e., if the trajectory is to 
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"be terminated by a condition of the form G(t f> x^)=:0, 
the program can handle it automatically.. 

Description 

DVDQ is called by: 

CALL DVDQCNEQ, T, Y, F, KD, EP, I FLAG, H, 

HMINA, HMAXA, DELT, TFINAL, MXSTEP, KSTEP , 

KEMAX, EMAX, KQ, YN, DT) 

Not all of the arguments must be supplied by the user. The 
arguments which must be supplied by the user will be listed and - 
briefly discussed below. Further explanation and definitions 
of the remaining terms may be found in the general comment 
section in the initial portion of the listing. 

NEQ = number of differential equations 
T = independent variable 
Y = dependent variable vector 

F = right-hand side of the vector differential equation 

KD = order of the differential equation (not to be 

confused with the order of o the integration scheme). 

For example, if xs f(t,x,x) is to be integrated 
directly, KD = 2 ; usually KD = 1 since x = f(t,x) 
is the usual trajectory analysis system. 

EP = local truncation error indicator. It is an absolute 
error indicator in the sense that the local error is 
kept less than EP in all components of the differential 
equation. If it is desired to control the error on 
each component of the vector differential equation 
separately, EP should be specified as a vector of 
negative numbers. The negative values alert the 
routine to a vector error specification. See the 
discussion of EP in the comment cards for further 
options . 

H c initial stepsise estimate. (Probably better to 

guess smaller than one would expect with a fourth- 
order scheme since the routine builds up from a 
first-order formula. However, the choice is not 
critical since the stepsize is adjusted rapidly.) 



HMINA = minimum allowable stepsize 
HMAXA = maximum allowable stepsize 

DELT ts output interval. Since the scheme only halves 
and doubles, use DELT = 2 n if possible; 
otherwise, it will interpolate for output values. 

TFINAL ss final value of the independent variable. 

MXSTEP = maximum number of steps allowable between 
output points. 

In addition DT(17,NEQ), YN( ) , Y( ) , F( ) , KQ( ) must appear 
in a DIMENSION statement in MAIN. See the example below or the 
comment portion of the listing for the particular dimensions. 
Operation 

Only one call is made* to DVDQ; thereafter the simple 
statement CALL DVDQ1 is used. The heart of operation of DVDQ 
is a computed GO TO statement which is driven by the parameter 
IFLAG, which has values 1,2, . . .,8. The full implications 
of each value of IFLAG are discussed in the initial comment 
section of the listing. Although the operation of DVDQ may 
appear complicated at first glance (because of the eight values 
for IFLAG) , the operation is straightforward with excellent error 
detection capabilities. Probably the easiest way to get 
acquainted with DVDQ is to study the simple program given below. 
Informal comments are given to explain roughly what each value of 
IFLAG is indicating. The listing follows the example. 

Example : Integrate x^=X2 

x^-x^+tx-j 

with t 0 =o,t f =2, x-j(0) = 1.0, x^ (0)-0.0. The following program 
will execute the integration with DVDQ. 
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t 

iistics* . 

"ISTICS* • 
' IN MAIN 


IMPLICIT REAL*8 (A-H,0-Z) 

SINCE SYSTEM IS FIRST_ORDE R WITH TWO DIFFERENTIAL EQUATIONS, 

Y f F, KQ, YN ARE OF DIMENSION 2. DT IS ALWAYS DIMENSIONED DT (17, NEQ) » 
DIMENSION Y <2) , F ( 2) ,DT (17,2) ,KQ (2) , YN {2) 

N EQ = 2 
T=0 . DO 
Y (1 ) = 1 e DO 
Y{2)=0.D0 
KD-1 

EP= 1 . D- 5 

H- 1 . D- 1 ' • • 

HMIN A = 1 . D- 4 
H MAX A=1 .DO 
DELT=2.D- 1 



. T.FIN AL = 2 . DO 
MXSTEP=1 000 

CALL DVDQ (NEQ,T ,Y ,F , KD , EP, I FLAG ,H , 3 MI NA , H MAXA , DELT ,TFIN AL , 
1MXSTEP, KSTEP, KEMAX, EM AX F KQ , YN , DT) 

GO TO 10 
CALL DVDQ 1 

GO TO (1, 1,3, 4,5, 6, 7, 8) , IF LAG 
F(1)=Y{2) 

F (2 ) =- Y (2 ) **2 + Y (1 -)*T 
GO TO 20 

WHITE (6 ,100) T, (Y(I) ,1=1,2) 

GO TO 20 

WRI TE (6 , 10 0) T, (Y(I) ,1=1,2) ^ 

STOP 0 
STOP 5 

STOP 6 , _ 

HMIN A=HMINA/2 , DO C, ^ cP Cawnot Be frrrA m &£> (PaoQ q£lV 

GO TO 20 _ Jf ^ Sm? o & 'buccc/tse' EP. 

STOP 8 


'XfLfiC = 1 <rx 2 - && 7B J)e*n/f)-riLfe' 
l-u & - n a*! <&~7^f+erj Call -^Vb^X, 
(£^Te.ub\t.TQ<z ) 2 — Cbe.de err c sS) 


3 ==£> Ou-r- PUT- Tb/^Tj 
C /QpTCtZ uj& mHGr y 

C =^> tf j OF 


a Jait&c fe^r/oftj Tj^/jse 7 - 
o*L St - of>< 

XFlP&zr r filXsTcP fl-rTfi>*ieh j Stop CjI XfJC&eAZU A\PSTEp. 


XdLnG ~7 =» tf<f 4 MtNA g&Quiceb r STOP oh £cserr /JMitJA oa.BP t 


FORMAT (* TIM E= * ,D23. 15/» Y (I) = ’ , 2D23. 15) 
END 


I LtU&ftL Pp£Atot,T£}l 

AfJ ti fiLLi hg S& Quebec Stop, 


SOURCE STATEMENTS = ’ 32, PROGRAM SIZE = 1400 

NO DIAGNOSTICS GENERATED 



SUBROUTINE DVDQ ( / N EQ/ , /T/ » / Y/ ,/F/»/KD/,/EP/ »/IFLAG/»/H/ , /HMINA/ 

* /HMAXA / » /D EL T/ ,/TF INAL/ ,/ MXSTEP/ , /KSTEP/ » /KEMAX/ , / EMAX/ , 

* /KQ/y/YN/ ,/DT/} 

DOUBLE PRECISION VARIABLE ORDER INTEGRATION SUBROUTINE 
FOR THE SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS. 

ANALYSIS AND CODING BY FRED T. KROGH , AT THE JET PROPULSION 
LABORATORY, PASADENA, CALIF.' APRIL 1, 1969. 

THIS DECK IN EBCDIC FORMAT . * 

CONVERTED FOR USE ON 360/75 BY MELBA W. NEAD, JPL APRIL, 1970. 

AT THE END OF THIS LISTING INSTRUCTIONS ARE GIVEN FOR REMOVING 
SOME FEATURES AND FOR ADDING OTHERS. THE GSTOP FEATURE IS ■ 
EXPLAINED NEAR THE END OF THE LISTING, i 

VARIABLES IN THE CALLING SEQUENCE HAVE THE FOLLOWING TYPES. 
INTEGER N6Q ,KD ( 1 ) , I F LAG ,MX STE P , KSTEP , KEMAX , KQ 1 1 > 

REAL EP (1 } ,HMINA, HMAXA, EMAX 1 

DOUBLE PRECISION T , Y ( 1 ) , F ( 1 ) * H , DEL T , TF I NAL , YN( 1 } ,D T [ 17 , 1 ) 

~ PARAMETERS WHICH- MUST BE ASSIGNED VALUES BEFORE CALLING 
DVDQ ARE NEQ, T, Y, KD, H, HMINA, HMAXA, DELT, 

TFINAL, AND MXSTEP. ! 

DVDQ IS USED ONLY ON THE INITIAL ENTRY. ALL OTHER 
ENTRIES ARE MADE BY CALLING DVOQl. IN ADDITION TO 
THE PARAMETERS MENTIONED ABOVE THE USER MUST ASSIGN 
•VALUES TO F (ONCE PER STEP INITIALLY, AND TWICE PER STEP 
AFTER GETTING STARTED) AND EP (EITHER INITIALLY, OR DURING 
THE INTEGRATION IF A RELATIVE ERROR TEST IS USED). 

THE FOLLOWING PARAMETERS GIVE ADDITIONAL INFORMATION ABOUT THE 
INTEGRATION AND ARE USED FOR STORAGE. THEY SHOULD NOT BE 
CHANGED BY THE USER. 1FLAG, KSTEP , KEMAX , EMAX, KQ,YN, AND DT. 


C AN EXAMPLE CF HOW ONE MIGHT SET UP THE CALLS TO DVDQ IS GIVEN 
C BELOW. ’ : 

IMPLICIT REALMS (A-H,0-Z) 

DIMENSION F (2) , DTI 17,2) ,Y(4) ,YNU) ,KQ(2) 

C SET PARAMETERS AND INITIAL CONDITIONS 

NEQ = 2 

KD = 2 ' . 

MXSTEP = 500 
EP = l.D-6 
HMINA = 0. 

HMAXA = 100. 

T - 0. 

Y( 1) = X. . 

Y { 2 ) = 0 • 

Y C 3 > = 0. 

Y(M = 1. 

H = 1. 

DELT = 1. 

TFINAL = 12. 

NEVALS - 0 



CC NOW HAKE FIRST ENTRY 

C CALL DVDQ(NEO,T»Y, F,KO,EP, IFLAG»H, HMINA,HMAXA, OELT, 

C • i TFINAL,MXSTEP,KSTEP,KENAX,EMAX,KQ,YN,DT) 

C GO TO 2 

C 

CC ALL SUBSEQUENT ENTRIES HADE HERE 

C 1 CALL DVDQl 

C 2 GO TO ( 10, 10,30,30, 50,60,50, 50) , IFLAG 
C '■ 

CC EVALUATE DIFFERENTIAL EQUATION , I FLAG = i OR 2 
C 10 R •= Y(l)**2 + Y(3)**2 ! 

C R = R*DSQRT(R) 

C F( 1) = -Y(l)/R 

C F( 2 ) = — Y (3 )/R 

C NEVALS = NEVALS+1 

C GO TO 1 

C 

CC OUTPUT, IFLAG « 3, OR FINISHED, IFLAG = 4 
C 30 PRINT 13, T,Y(X),Y(3),KSTEP, NEVALS 
C 13 FORMAT ( * 1 , F6 *2 , 1 P2 D2Q . 12, 2 1 5 i 
C IF 1 IFLAG. EQ. 3) GO TO 1 

C STOP ■ 

C 

CC ERROR CONDTION, IFLAG = 5, 7 OR 8. 

C 50 PRINT .12, 1 FLAG, H .... 

C 12 FORMAT ( * I F LAG, , 13 , D20 . 10 ) 

C GO TO 30 

c 

CC , EP TOO SMALL, IFLAG = 6 

C 60 EP. = -32.*EMAX*EP ‘ 

C PR INT 14, T , E P - 

C 14 FORMAT { 5 T , NEW £P= * ,F6 .2,015.6) 

C • GO TO 1 ■ / - . • ' 

C END 

C 

C 

C • « *•*•«•••••*«*• ........ «>••••*»••»•«»•••<>••••••••« »«•«»#«# » . . . « 

c • .■ ■ . . .... . 

C ■ ■ . ■ 

C THE USAGE OF THE VARIABLES IS GIVEN BELOW. 

C 

C NEQ=NUMB ER OF EQUATIONS (INPUT) 

C - ... 

C T = INDEPENDENT VARIABLE (INITIAL VALUE SUPPLIED BY THE USER) 

C 

C ‘ Y = CURR ENT VALUE OF DEPENDENT VARIABLE. THE INITIAL 

C VALUE OF Y MUST BE SPECIFIED BY THE USER BEFORE 

C-- THE. FIRST ENTRY. ...THE . DIMENS ION OF Y MUST BE 

C AT LEAST AS GREAT AS THE SUM OF THE ORDERS OF 

C THE DIFFERENTIAL EQUATIONS WHICH ARE BEING 

C INTEGRATED. IF WE LET KD(I) DENOTE THE ORDER 

C OF THE I-TH DIFFERENTIAL EQUATION, THEN Y(J) 

C • IS THE K-TH DERIVATIVE OF THE L-TH COMPONENT , 

C WHERE L IS THE SMALLEST INTEGER FOR WHICH 

C KD( 1) + KD ( 2) +. *.+KD(L) • GE • J AND K-KD { L )+ J-l- { KD( 1 ) 

C +KD(2 > + . . .+KD(L ) ) , J = I, 2 ,. . . »( KD{ 1 ) + KD{ 2) + . . . + KD(NEQ) ) . 

C (FOR EXAMPLE, FOR THE SYSTEM F(1)=UPP, F(2)=VPP, WHERE P 

C DENOTES A PRIME, Y(1)^U, Y(2)=UP, Y(3)=V, Y(4)=VP.) 

C 

C F ( I ) =KD ( I ) — TH DERIVATIVE OF THE I-TH COMPONENT WITH RESPECT 



c ■ TO T, I= 1 , 2 *...»NEQ. THE USER MUST PROVIDE 

t THE CODE WHICH COMPUTES F GIVEN Y AND T. 

C 

C KD GIVES THE ORDER OF THE DIFFERENTIAL EQUATIONS IN THE 

C SYSTEM. KD MUST BE LESS THAN OR EQUAL TO 4. 

C (FOR DIFFERENTIAL EQUATIONS WITH DIFFERENT ORDERS SET 

C KD.LT-0 « IF THIS IS DONE IT IS ASSUMED THAT KD IS A VECTOR 

C AND THAT A8S(KD(I)) GIVES THE ORDER OF THE I-TH EQUATION.) 

C 

C EP IS A PARAMETER USED TO CONTROL THE LOCAL ERROR. 

C IF EP IS POSITIVE THE LOCAL ERROR IS KEPT LESS 

C THAN EP IN ALL COMPONENTS OF THE DIFF* EQ. 

C (THE ESTIMATED LOCAL ERROR IS KEPT LESS THAN EP IN 

C THE- ( KOI I) —1 ) -ST DERIVATIVE OF THE I-TH . COM PON ENT • THUS 

C FOR EQUATIONS WITH ORDER GREATER THAN ONE, THE ERROR 

C IN A DERIVATIVE IS ESTIMATED. IN THIS CASE THE VALUE OF 

C EP REQUIRED TO OBTAIN A GIVEN ACCURACY IN THE DEPENDENT 

C VARIABLE DEPENDS ON THE SCALING.) 

C IF EP.LT.O, THEN IT IS ASSUMED THAT EP 

C ISA VECTOR. LET K BE THE SMALLEST VALUE 

C OF I FOR KH ICH EPIIKGE.O. FOR I.LT.K 

C THE LOCAL ERROR CONTROL IS BASED ON . 

C ’ A8 S ( EP ( I ) ) » AND FOR I.GE.K IT IS BASED ON 

C EP(K). IF ONE WANTS A RELATIVE ERROR TEST — ; 

C • FOR EXAMPLE, THE LOCAL ERROR IS TO BE KEPT 
C LESS THAN C*P WHERE C IS A CONSTANT 

C AND P IS A POSITIVE FUNCTION OF T AND Y, 

C THEN ONE SHOULD SET EP=C*P WHEN IFLAG=1. 

C. ‘ IF EP = Q AND HMAXA.NE.O, I FLAG IS SET EQUAL 8. IF EP = 0 AND HMA XA =0 , 

C NO ERROR TESTS ARE PERFORMED AND THE ORDER IS ) AND STEPSIZE ARE 

'C ‘ NOT CHANGED. THIS OPTION SHOULD NOT BE USED IF KQ(I)=L FOR ANY I • 

C 

C IFLAG IS USED FOR COMMUNICATION BETWEEN THE INTEGRATOR .. ....... 

C ........ AND THE PROGRAM WHICH CALLS- IT . THE VALUE 

C OF IFLAG SHOULD NOT BE CHANGED BY THE USER^. 

C THE FOLLOWING VALUES OF IFLAG HAVE THE FOLLOWING MEANINGS. 

-C - =1 -THE ■ VALUE -OF- -Y -FOR THE CURRENT STEP HAS BEEN • . 

C - PREDICTED. THE USER SHOULD COMPUTE F AND CALL DVDQ1 ^ 

C IF A RELATIVE ERROR TEST IS USED THE NEW VALUE 

C OF EP SHOULD ALSO BE COMPUTED HERE. 

C =2 THE VALUE QF. Y FOR THE CURRENT STEP HAS BEEN 

C ■ CORRECTED. THE USER SHOULD COMPUTE F AND CALL DVDQX . 

C =3 AN OUTPUT POINT HAS BEEN REACHED (SEE DESCRIPTION 

C OF DELT ) t PRINT RESULTS AND CALL DVDQ1 . 

C =4 T = TF I NA L IF DVDQI IS CALLED WITH T^TFINAL AND 

C I FL AG=4 * IFLAG IS SET EQUAL TO 8. IF THE VALUE OF 

C TF INAL IS CHANGED THE INTEGRATION WILL CONTINUE. 

C =5 KS TE P= K SO UT (SEE THE DESCRIPTION OF MXSTEP). 

C =6 EMAX.GT..1 AND IT APPEARS TO THE SUBROUTINE THAT 

C REDUCING H WILL NOT HELP BECAUSE OF ROUND- OFF ERROR. 

C IF THIS OCCURS A LARGER VALUE OF EP (OR OF AB S ( E P I KEM AX > ) IF 

C EP IS A VECTOR) SHOULD PROBABLY BE USED. IF EP IS NOT 

-C INCREASED, TOO SMALL A STEPSIZE IS LIABLE TO BE USED. (WE HAVE 

C FOUND THAT REPLACING EP WITH 32.*EMAX*EP WORKS QUITE WELL.) 

C INCREASING EP IN THIS WAY- WILL NOT DEGRADE THE ACCURACY, 

G- HOWEVER IF. THE NATURE OF. .THE- PROBLEM CHANGES IT MAY PAY TO 

C USE A SMALLER VALUE OF EP LATER IN’ THE INTEGRATION. 

C =7 A8S (H) . LT.HMINA. TO CONTINUE WITH THE CURRENT 

C VALUE OF H, SET HM I N A. L E . ABS ( H ) AND CALL DVDQI. 

C IF THE INTEGRATOR HAS JUST HALVED H ONE MAY CONTINUE 













WITH TWICE THE STEPSIZE BY SIMPLY CALLING DVDQl * (SUCH 
AN ACTION IS RISKY WITHOUT A CAREFUL ANALYSIS OF THE 
SITUATION,) IF THE STEPSIZE HAS NOT JUST BEEN HALVED 
(ABS(H) .LT.HMIN A MAY BE DUE TO THE USER INCREASING THE 
VALUE OF HMINA OR TO HAVING TOO SMALL AN H AT THE END 
OF THE STARTING PHASE.) THE INTEGRATION WILL CONTINUE 
WITH THE CURRENT VALUE OF H AND A RETURN TO THE USER WITH 
TFLAG^ WILL BE MADE ON EVERY STEP UNTIL ABS(H)*GE, HMINA. 

=8 ILLEGAL PARAMETER IN THE CALLING SEQUENCE. IF DVDQ1 
IS CALLED WITH I FL AG= 8 THE PROGRAM IS STOPPED. 

H=CURRENT VALUE OF THE STEPSIZE. IN SELECTING THE INITIAL 
VALUE FOR H, THE USER SHOULD REMEMBER THE FOLLOWING — 

1. THE INTEGRATOR IS CAPABLE OF CHANGING H QUITE QUICKLY AND 
THUS THE INITIAL CHOICE IS NOT CRITICAL. 

2. IF IT DOES NOT LEAD TO PROBLEMS IN COMPUTING THE DERIVATIVES 
(E.G. BECAUSE DF OVERFLOW OR TRYING TO EXTRACT THE SQUARE 
ROOT OF A NEGATIVE NUMBER } » IT IS BETTER TO CHOOSE H MUCH 

• • TOO LARGE THAN MUCH TOO SMALL. 

3. IF H*DELT.L£.0 INITIALLY, AN IMMEDIATE RETURN IS MADE 
WITH I F LAG=8 • THE SIGN OF H IS WHAT DETERMINES THE 
DIRECTION OF INTEGRATION. 

, . k, IF DELT=H*C2**K) K A NONNEGATIVE INTEGER THEN OUTPUT 

. VALUES WILL BE OBTAINED WITHOUT DOING AN INTERPOLATION. 

HMINA AFTER GETTING STARTED, AND WHENEVER H 

IS HALVED* ABS { H ) IS COMPARED WITH HMINA. 

IF AB5(H) 0 LT. HMINA CONTROL IS RETURNED TO 
THE USER WITH IFLAG-7. 

HM AX A THE STEPSIZE IS NOT DOUBLED IF 

DOING SO WOULD MAKE ABS l H ) .GT .HMAX A 

DELT ENABLES THE USER TO SPECIFY THE POINTS WHERE 

OUTPUT IS DESIRED. LET TOUT -DELT + THE VALUE OF T THE LAST 
TIME CGNTRGL WAS RETURNED TO THE USER WITH I F LAG=3 • (INITIALLY 
TOUT =T HE INITIAL VALUE OF T.) CONTROL IS RETURNED TO THE 

.-USER WITH IFLAG=3. WHENEVER T-TOUT..IF TOUT DOES NOT FALL . 

ON AN INTEGRATION STEP, OUTPUT VALUES ARE OBTAINED BY 
INTERPOLATION ON THE FIRST STEP THAT ( T-TOUT ) *H. G T . 0. 
INTERPOLATED VALUES FOR BOTH Y AND F ARE COMPUTED. 

(NOTE THAT A RETURN WITH IFLAG=3 IS ALWAYS MADE 

BEFORE TAKING THE FIRST STEP.) ... ' ; 

TFINAL CONTROL IS RETURNED TO. THE USER WITH I FLAG=4 WHEN 
T REACHES TFINAL. IF TFINAL DDES NOT FALL ON AN INTEGRATION 
STEP VALUES AT TFINAL ARE OBTAINED BY EXTRAPOLATION. 

MXSTEP ON THE INITIAL ENTRY, AND ON ENTRIES 

WITH 2.LT. IFLAG.LT.6 KSOUT IS SET EQUAL TO 
KSTEP+MXSTEP. AT THE END OF EACH STEP KSTEP IS INCREMENTED 
AND COMPARED WITH KSOUT. IF KSTE P. GE . K SOUT CONTROL IS 
RETURNED TO THE USER WITH IFLAG^S. {THUS IF DELT IS 
SUFFICIENTLY LARGE, CONTROL WILL BE RETURNED TO THE USER 
WITH I F L A G= 5 EVERY MXSTEP STEPS.) 

KSTEP -NUMBER OF INTEGRATION STEPS TAKEN (COMPUTED. , . .- i; . 

BY THE . INTEGRATOR. ) . .. ‘ ' - 

KE MAX = I NO EX OF COMPONENT RESPONSIBLE FOR THE 


on. . ■ ' o o r> nonnnnnnnnn o n n n n n no ono 


C VALUE OF EMAX (SEE BELOW) . 

•c 

C EM AX=L ARGES T VALUE IN AMY COMPONENT OF (ESTIMATED ERRORJ/EP .. . 

C- -ORDINARILY- THE.. -STEPS I ZE IS -HALVED IF. EMAX .GT*. 1* WITH A 

C RECENT HISTORY OF LOCAL ROUND-OFF PROBLEMS VALUES OF EMAX AS 

C LARGE AS I ARE PERMITTED. THE STEPSIZE IS NOT HALVED ON ANY 

C STEP THAT ROUND OFF ERROR APPEARS TO BE LIMITING THE PRECISION. 

C 

C KQ ( I ) = H I GHE ST ORDER DIFFERENCE USED IN INTEGRATING 

C THE I -TH EQUATION. [COMPUTED BY THE INTEGRATOR) 

C 

C YN=A VECTOR WITH THE DIMENSION OF Y USED TO STORE 

. C THE VALUE OF Y AT THE END OF EACH INTEGRATION STEP. 

C 

C DT =AN ARRAY WITH DIMENSION DT { 17 » NEQ ) USED TO 

. C STORE THE DIFFERENCE TABLE. 

c • 

c < 

DOUBLE PRECISION TO UT S TL , TPD ,TPD1 , TPD2 ,HH , F AC 
DOUBLE PRECISION DO , D, GAM , GA S, PT , TP 

DIMENSION DO (19) tD { 18) ,GAM<17*4) ,GAS ( 17 ) t PT ( 18 ) ,FAC(3> 

EQUIVALENCE ( DD ( 2 ) * D ( 1) ) 

' DIMENSION E TA (1 5 * 15 ) . 

C ' - - - ■ - 

DATA KMAXO/4/ 

KM A XO IS THE MAXIMUM ORDER DIFFERENTIAL EQUATION 
THIS SUBROUTINE WILL INTEGRATE. 

DATA FAC/1. DO» . 5D0 , . 166666666666666666700/ 

FAC t J J = 1/ (FACTORIAL J)t J= l, 2, . . . ,MAX( 2 , KMA XO- 1 ) 

DATA KQMAX/ 16/ ' '- 

KQMAX GIVES THE MAXIMUM ORDER. - 

THERE IS LITTLE POINT IN HAVING KQMAX MUCH BIGGER THAN THE NUMBER 
OF DECIMAL DIGITS IN THE MANTISSA. * 

IF KQMAX IS SET LESS THAN 6 t DT, 0, AND PT SHOULD BE DIMENSIONED 
AS IF KQMAX-6 . . . 

DATA R ND»KBIT2/8. 88 E— 16,108/ 

RND IS APPROXIMATELY 2**<3-B) WHERE B IS 
THE NUMBER OF BITS IN THE MANTISSA. 

KB IT2= 2*8+2 WHERE B IS THE NUMBER OF BITS IN THE MANTISSA. 

IF THE DERIVATIVES ARE NOT COMPUTED TO THE ACCURACY EXPECTED 
FROM THE WORD LENGTH OF THE COMPUTER (FOR EXAMPLE BECAUSE OF 
. CANCELLATION PROBLEMS OR TABULAR DATA), THEN THESE CONSTANTS 
CAN BE CHANGED TO REFLECT THE NUMBER OF BITS WHICH ARE 
SIGNIFICANT IN THE COMPUTED DERIVATIVES. (THIS IS NOT NECESSARY, 
BUT IS WISE IF THE ACCURACY REQUESTED IS DIFFICULT TO OBTAIN 
BECAUSE THE DERIVATIVES HAVE SO FEW SIGNIFICANT DIGITS.) 

DATA PI » P01 , P25 , P3E1/ • 1, *01, .25,3./ 

THE ABOVE DATA STATEMENT CONTAINS VARIOUS CONSTANTS 
USED IN THE SUBROUTINE. 

.. DATA PT/ 1. DO, 2-DO, 4,. DO, 8. DO.,16. DO, 32.00,64. DO, 128. DO ,256. DO , . . .. 

1 512. 00,10 24 .DO, 20 48.00,4 096.00,8 19 2. DO, 16 384. DO, 32768. DO, 

2 65536.00,131072.00/ 

PT{ J) = 2**< J-l) , J = 1 ,2,...,KQMAX+2 

DATA PTS1, PTS2,PTS3 ,P T S4 ,PTS 5 ,P 5/ 1 . ,2. ,4. ,8. , 16. 5/ 

DATA PT S 1 , P TS2 , PTS4,PTS5 ,P5/1. ,2. , 8. ,16. ,.5/ 


f- V* 



o o o 


BIO 


C 

DATA GAS/1 . 00,-*5D0 »-8 . 33333 3333 333 3 33 3 3 3D- 02, 

1 -A. 1 6666666666666667D-02 ,-2»63888888883888809D-02, 

2 —1 <,0750-0 2? — 1.42 69179894179894 20— 02, 

3 -1.136739417989417990-02,-9.356536596119929450-03, 

4 -7. 892554 01 234567-901 D-03 , -6. 7858499846347068 6D-03 , 

5 -5 .924056 4123376 6 2340-03, -5. 236693 2579502850 7D-03, 

6 -4. 677498407042264520-03 ,-4.214952239005472860-03, 

7 -3.82 6 89955321 188 442D-03 ,^3. 497349 84534991 76 5D-03/ 

G AS { I > GIVES THE I-TH ADAMS-MOULTON CORRECTOR 
COEFFICIENT, 1 = 1 , 2 , . . . , KQMA X+ 1 . • ■ • • 

DATA GAM! 01, 01) , GAM (02 ,01 ) , GAM ( 03,01 ) , GAM { 04,01 ) , GAM { 05 ,01 ) , 

* GAM { 06 ,01), GAM! 07, 01) ,GAM{ 08,01) ,GAM( 09,01) ,GAM ( 10 ,01) , 

* GAM (11 ,01) , GAM {,12,01 ) , GAM ( 13,01) , GAM ( 14,01) , GAM ( 15,01), 

* GAM( 16,01) , GAM (17 ,01)/ 

1 1. DO, „5D0,.4166666666666666667D0, .37500, 

2 .34861111 1111 1111 11 DO,. 32986111111111 1111D0, 

-3 ,-.31559193 12 16931 2 17D 0 , . 3 042245 3 703703 7 037D0 

4 *29486800 04 40 917 108D0, .2869 7544642 857 1429D0,. 

5 : .28018959644393672200, . 2 742655 400 3 1599 059D0 , 

-6 .269028846773648774D0, .2643513483666 06.51000, 

7 -.26013639612760103700, . 2563094965743 89153D0, 

8 . 2528 1214672903 92 350 0/ 

DATA GAM (01,02) , GAM (0 2,0 2) ,GAM( 03,02) P G AM l 04,02) ,GAM(05 ,02) , 

* GAM( 06 ,02 ) ,GAM (07 ,02) , GAM ( 0 8 ,02 ) , GAM ( 09, 0 2 ) , GAM ( 10,02), 

* GAM ( 11 ,02) ,GAM( 12 ,0 2) , G AM { 13,02) , GAM (14 ,02) ,GAM( 15 ,02 ) , 

* GAM{ 16 ,02 ) ,GAM ( 17 ,02)/ 

1 . 5D0, . 166 666666 66 66666667D0 , .125 DO,. 10 555555555 5 5555556D0, 

2 9 .3750-2 , 8. 5615079365 07 93651D -2, 

3 7. 957 1759 2592592 5 93D-2,7 .485 22 927689 594356D-2 r 

4 7. 1 03 2986 11111111 110-2,6. 7858 4998 463470686 D-2 , 

5 6. 516 46 20535 7 142857 D-2, 6. 2840319095403420 8D-2, ' 

6 6. 0807479 291 54943 87D-2 ,5 . 9 009331 3460766200 D-2 , 

7 5.740349329817026630-2,5. 5957 5975 2 559 8682 5D-2 , 

B 5.464643933250064670-2/ 

DATA GAM( 01 ,03) ,GAM{ 02 ,03) , GAM (03, 03 ) ,-GAM( 04 ,03 ) , GAM ( 05 ,03 ) , 

* GAM (06, 03) , GAM ( 07,03) , GAM ( 0 8,03 ) ,GAM{ 09 , 03 ) ,G AM ( 10 ,03 ) , 

* GAM{ 11,03) , GAM (12 ,03 ) , GAM { 13,03 ) ,GAM[ 14,03) , GAM t 15,03) , 

* GAM( 16,03) ,GAM(17, 03) / 

1 .1666666666666666667D0, 4 .166666666666666670- 2, 

2 2.91666666666666667D-2,2.36111111111111111D-2, 

3 2.0337301 5 873 015 873D-2, 1.8 1299603 1746031 750- 2, 

4 1*651813271604938 27D-2 , 1.52772266313932981D-2, 

5 1.4285188 19143819 14D- 2, 1 . 3469 396 5 53 1639143D-2 , 

6 1. 2783 6579217 0975 70D-2,1 . 21970388 2312 38926D- 2, 

7 1. 168796164557332160-2,1. 1 24086 63 3 52 884755 D-2 , 

8 1.08442182943468791D-2, 1.048926554478428630-2, 

9 1.01692338611494262D-2/ 

DATA GAM( 01,04) ,G AM (02,0 4) , G AM ( 03,04-) , G AM ( 04 , 04 ) , G AM ( 05 ,0 4 ) , 

* GAM (06 ,04) , GAM (07,04) , GAM ( 0 8 ,04) ♦ GAM( 09,04) ,GAM( 10,04) , 

* GAM( 11,04) ,GAM( 12 ,04) ,GAM( 13 , 04 } , GAM ( 14 ,04) , GAM ( 15 ,04) , 

* GAMl 16, 04) ,GAM( 17,04)/ 

1 4. 166666666666666670-2, 8. 3 333 3 33 3 333 3333 3 3 D-3 , 

2 5. 55 55 55 5 55 55555 5 56D-3, 4. 365 0793 650 7 *5365 08D-3 , 

3 3.6789021 16402 116 400- 3, 3 . 2 2365 5202 82 186949D- 3 , 

4 2.895447530864197530-3,2.645435339880284320-3, 

5 2.4473749 148 2 283 1 490-3, 2 . 2 857 5543 8 1 80 5241 60- 3 , 

6 2. 1509366 94 81 4836 35D-3, 2. 036308 710202283840-3, 

7 1.9374130 1 123 4333 02 D-3 , 1 . 85102419 106078320D-3 , 
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B 1 .7747637 478 L 296400 D~3, X .706835646 05258723D- 3 ? 

9 1.6458559 10054651 58D-3/ 

GAM(ItJ) GIVES THE I-TH ADAM S- FALKNER PREDICTOR 
COEFFICIENT FOR INTEGRATING J-TH ORDER DIFFERENTIAL 
EQUATIONS p 1=1 ,2 , • » KQMAX+- 1 * J"1 ,2 » - . - pKMAXO. 


DATA 


1 . 
2 . 

1. 
9. 

DATA 


3. 

8 . 

4. 
2 . 

DATA 


3 

U 
8 . 
5. 

DATA 


* 

❖ 

1 

2 

3 

4 

* 

* 


2 . 
2 . 

1. 
8 . 
DATA 


ETAt 
ETA( 
ETAt 

1363636 OE 
713 815 9 OE 
40603000E 
031 3 726 OE 
ETAt 
E TA { 
ET A { 

40909090E 
1414478 OE 
218 0901 OE 
709 4118 OE 
ET A ( 
ETAt 
ETAt 

42857140E 
503 0365 OE 
0344571 OE 
23 1 9680 OE 
ETAt 
ETAt 
ETAt 

85 71428 OE 
0875506 OE 
205 1 685 OE 
0970932 OE 
ETAt 
ETAt 
ETAt 


1 2 . 42424240E 

2 2 . 2773280 OE 

3 1.5223182 OE 

4 1.0796 1240E 

DATA ETAt 

* ETAt 

ETAt 


2* 097902 10E 
1.912 9555 OE 
1 .657 63530E 
1 • 2789254 OE 
DATA ETAt 


* 
* • 
* 
1 
2 


ETAt 

ETAt 

1. 84615380E 
6.961 9269 OE 


01 , 01 ) 

06,01) 

11 , 01 ) 

- 01 , 6 
-02, 2 
-02, 1 
-03/ 

01 , 02 ) 
06,02) 
11 , 02 ) 

- 01 , 2 
- 02 , 6 
-02, 3 
- 02 / 

01.03) 

06.03) 

11.03) 

- 01 , 3 
-01, 1 
-02, 7 
- 02 / 

01 ,04) 

06.04) 

11.04) 

- 01 , 2 
-01, 1 
01 , 1 
-02/ 

01 .05) 

06.05) 

11.05) 

-01, 2 
-01, 2 
- 01 , 1 
- 01 / 

01.06) 
06,06) 
11 ,06) 

-01, 2 
- 01 , 1 
- 01 , 1 
-01/ 

01.07) 

06.07) 

11.07) 

-01, 2 
-02, 1 


,ETA(02,01),ETAt03,01),ETAt 04 ,01), ETAt 05, 01) , 
,ET At 07 ,01 ) ,ETA( 08,01 ), ETAt 09,01 ), ETAt 10 ,01 ) , 

, ETA t 12,01) ,ETAl 13,01) ,ETAt 14,01) ,ET At 15,0 1 )/ 
3. 33333 33 OE— 01, 2 . 5 OOOOOOOE-Ol , 

. 73076930E-02 , 4. 60526330E-02 , 3 .43749980E-02 , 
.225473 10E-02, 1 . 87484580E- 02 , 1.611232206-02 , 
*241970 60E— 02 1- 10802170E-02 , 9 ;96793590E-03, 

, ETA (02, 02 ), ETAt 03,02) , ETAt 04,02) , ETAt 05 ,02) , 

,E T A ( 07 ,02) , ETAt 08,02 ) , ET A ( 09 ,02 ) , ETA t 10 ,02 ) , 
,£TA{ 12,02) , ETAt 13,02) » ETA t 14,02) , ETAt 15 ,02) / 
2.00000000E-01, 4.00 00 0000 E- 01, 

. 01923080 E- 01 , 1 . 38 1 57 900E-01 , I .03124990E-0 1 , 
.676419 30E-02 , 5. 62453730E-02, 4. 833 6967 OE- 02 , 
o 72 5 91 1 70 E— 02 , 3 . 324 0651 OE-02 , 2 . 9903808 0E-02 , 

, ETA (02 ,03) , ETA ( 03, 03 ), ETA (04,03) , ETAt 05,03), 

, ETA (07,03) , ETA tO 8,03) , ETAt 09 ,03), ETA (10,03 ) , 

, ETA (12 ,03), ETA (13, 03), ETAt 14,03) , ETA £15,03) / 
I. 42 857140 £-01 , 2 .85714280 E-0 1 , . 

♦ 46153840E— 01 , 2 . 45614040E- 01 , 1 . 875 OOOOOE-Ol t . 
.24626490E-01 , 1 . 058 73640E- 01 , 9 . 15 858320E- 02 , 

• 12783130E — 02 , 6 . 38220510E-02 , 5. 75925 170E-02 , 

, E T A ( 02 ,04) ,ETAt 0 3,04) ,ETA(04,04) , ETAt 05 ,04) , 

, ETA ( 07, 04), ETA (08, 04) ,ETA(09,04),ETA(I0,04) , 

, ET A 1 12,04 ), ETA ( 13, 04), ETAt 14,04) ,ET At 15,04)/ 
l.limilOE-Ol, 2*2 2222220E-0 1 , 
.53968250E-01 , 3 . 07C17 540E-0 1 , 2. 50000000E-01 , 
.780 378 50 E-0 1 , 1 . 54399060E-0 1 , 1 . 3 5 68 271 GE-0 1 , 
.07997450E-01 , 9. 75 05 9080E-G2 , 8 . 86038 72 0E-02 , 

,ETAt02 ,05) , ETAt 03,05 ) , ETAt 04,05 ) , ETAt 05,0 5 ) , 

, ETA (07,0 5) ,ETA( 0 8,05) ,ETA( 09,05) ,ETAQO ,05) , 

, ETA (12., 05 ) , ETAt 13,0 5 ), ETAt 14,05) , ETAt 15,05)/ 
9. 0909091 0E-02 , 1. 8 18 18 180 E-0 1 , 
.42424240E-01 , 1 . 73 1601 70E- 01 , 2. 50000000E-01 , 
.054282 90 E— 01 , 1. 85278880 E-0 1 , 1 .67608050E-01 , 
.38853860E-01, 1 . 27 1 81 62 0E-01 , 1 . 16957 100E-0 1 , 

,ETA(02 ,06) ,ETA 103,06) ,ETAt 04,06 ) , ETA (05 ,06) , 

, ETAt 07,06) , ETAt 0 8,06) ,ETA( 09,06) , ETAt 10 ,06) , 
,ET A t 12,0 6) , ETAt 13, 06), ETAt 14,06) , ETAt 15,06) / 
7. 69230 76 OE- 02 , 1 . 53 846 150E-0 1 , 

• 23 7762 20 E-0 1, 1 . 86480190E- 0 1 , 1. 1 188811 0E-01 , 

.91733070E— 01 , 1 . 852 78880E-01 , 1 . 7598 8460E-0 1 , 

* 555 163 30E-0 1 , 1 . 456 80770E- 01 , 1 . 36449950E-01 , 

,ETA( 02,07) ,ETAl 03,07 ) ,ET At 04,07) , ETA (05, 07) , 
,ETA{ 07,0 7) , ETAt 08,07) , ETAt 09,07) ,ETA( 10,07) , 
,ET A (12 ,07) , ETAt 13,07) , ETAt 14,0 7) , ETAt 15,07)/ 
6. 6 6 6666 6 OE- 02 , 1 . 3333 33 30E -01 , 

.05 128200 E- 01 , 1 . 86480190E-0 1 , 1. 34265730E-01 , 

.3 94 42230 E— 01 , 1 . 52023690E-01 , 1. 56434 190E-0 1 , 
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3 1.56012730 E— 0 I * 1 .5270 79 70E-0 1 , 1 . 4799 3160E-01 , 1 . 423 82 55 OE-O 1 , 

4 1.3641 372 0E -01/ 

DATA ETA( 01,08 > , ETA { 02 , 08 ) , E TA ( 03 , 08 ) , E TA ( 04 , 00 ) ,E T A { 05 ,08 ) , 

* ETA{ 06 ,03 1 ,ETA [07 ,0 8 5 ,ETA( 08,08 ) ,ETA(09 f OB) , ETA ( 10 ,08) , 

* ETA ( 11,08) ,ETA( 12 ,08) ,ETA{ 13,08) ,6TA( 14,08 ) ? ETA (15 ,08 )/ 

* 5. 88235290 E- 02, 1 . 1 76470606-01 , 

1 1. 6470588 OE -01, 1 . 88 2352 90 E-01 , 1 . 80995470E-01 , 1 .44796380E-0 1 , 

2 9.21431500E-02, 4 . 2 i 2 258 30E- 02 , 9 . 772 551 90E- 02 , 1 . 14931240E-01 , 

3 1.25367360E-01, l .30961120 E-01 , 1 . 33 1 93 8 50 E-01 , 1.33136920E-01, 

4 1. 31546610E-0I/ 

DATA ET-A(01, 09 ) ,ETA(02, 09), ETA ( 03,09) ,ETA( 04,0.9) »E.TA( 05,09) , - - 

* ETA (06,09) , ETA (07 ,09) , ET A ( 0 8 , 09 ) , ET A ( 09 , 09 ) , ET A ( 10 , 09 ) , 

• * ETA( 11 ,09) ,ETA( 12,09) , E TA ( 13,09) ,ETA( 14,09) ,ETA( 15 ,09) / 

* • 5. 26315790 E- 02, 1 .0 5263 160E-0 1 , 

1 1. 486068 10E— 01, 1 .733746 1 OE-Ol , 1 . 73 3746 10E-01 , 1 *48606 81 OE-O 1 , 

2 1.06692070E-0X , 6 .096 689 70E- 02 , 2 . 4941 00 30E-02, 6. 63064 84 OE- 02 , 

3 8. 3578253 0E-02, 9 .62 9493 90 E-02 , 1 . 05 1530 40 E-01 , 1 . 10947 43 OE-Ol , 

4 1. 1438835 OE-Ol/ 

DATA ETA ( 01, 10 ), ETA (02 ,10) , ETA ( 0 3, 10),ETA( 04 , 1 0) , ETA ( 0 5 , 10 ) , 

* . ETA(06,10) ,ETA(07,10) ,ETA(08,10) ,ETA( 09 ,10 ) , ETA l 10,10) , 

* ETA( 11,10) .,ETA( 12,10) ,ETA{ 13,10) ,ETA(14, 10) ,ETA( 15 ,10) / 

. * 4.76190480E-02, 9 .52380950E-02, 

1 1.3533 83 5 OE— OX , 1 . 60401 OOOE-Ol , 1. 65 1 18680E-01 , 1 *4860681 OE-Ol , 

2 1.15583080 E— 01, 7 .548 28 240E-02 , 3 . 9 193 00 50E- 02 , 1 . 451 5928 0E-02 , 

.3 4.37790860E-02, 5 . 88469 080 E-02 , 7 . 1400209 0 E-02,.. 8 . 136 1 434GE- 02 » . . . 

4 8 .896 8707 OE-O 2/ 

■ DATA ETA (01 ,11) , ETA (02, 11) , ETA{ 03, 11 ) , ETA ( 04,11 ) ,ETA ( 05 ,11) , 

ETA(06,11) ,ETA(07,11) , ETA ( 08 , 1 1 } , ET A ( 09 , 1 1 ) , ETA { 10,11), 

ET A( 11,11) ,ETA( 12,11) ,ETA( 13,1 1) ,ETAt 14,11) ,ETA( 15,11) / 

4.34782610E-02, 8 . 6 956 5 2 10 E-0 2 , 

1 1 • 2422360 OE-Ol * 1 .49 06 8 320E-0 1 , 1 . 5691402 OE- 01 , 1 . 46 45309 OE-O 1 , 

2 1 .206 0843 OE-O 1 , 8 .61 4887 60 E-02 , 5 . 16 893250E- 02 , 2 . 46 1 3 964 OE-O 2 , 

3 8.3308803 0E-03, 2 . 824651 60 E-02 , 4 . 03201 I80E-02 , 5 . 1386 1 540E-0 2 , 

4 6.I0071000E-02/ 

DATA ETA (01,12) ,ETA (02 ,12) , ETA ( 03 , 12 ) , ETA( 04 , 12) , ETA (05,12) , 

* ET A ( 06, 12)., ETA ( 07, 12) , E TA ( 0 8 , 1 2 ) , E TA ( 09 , 1 2 ) , E T A ( 10 , 1 2 ) , 

* - • ETA (1 1 ,12 ) , ETAX 12,12) ETA ( 13 , 12 ) , ET A ( 14,12) , ETA t 15,12) / . 

* ■ .. . 4.00000000E-02 , 7 . 9 9999990 E-0 2 , 

1 - 1 . 1478 261 OE-Ol , 1. 39 1 30430 E- 0 1 , 1 . 4906 8320E- 0 1 , 1- 43105590E-01 , 

2 1.2302059 OE-Ol , 9 .3 72997 70 E-02 , 6 . 2 02 71900E-02 , 3 . 44595 500 E-02 , 

3 1.51622020E-02, 4 .725 881 20E-03 , 1 . 78691420E- 02 , 2 . 699 06 87 0E-02 , 

4 3.60496320E-02/ 

DATA ETA ( 01 ,13) ,ETA( 02 ,13 ) , E TA t 0 3 ,13) , ET A( 04 ,13 ) , ET A ( 05 ,13 ) , 

* ET A (06, 13) , ETA (07 ,13) , ETA ( 08,13) , ETA ( 09,13) , ETA ( 10,13) , 

* ETA 111, 13) ,ET A { 12 ,13) , ETA ( 13,13) , ETA (14,13 ) , ETA t 15,13)/ 

* 3. 7037037 OE- 02, 7. 4 0740740E-02 , 

1 1.06666670E-01, 1 .30370370E-01, 1 . 41 70692 OE-O 1 , 1. 39 130 430E-0 l , 

2 1. 2367150 OE— 01, 9 . 8 93 71 9 80E-02 , 7 . 029 74820E-02 , 4 . 3393508 0E-02 , 

3 2. 2462522 0E-02, 9.18921340E-03, 2 . 65 466170E- 03 , 1. 11137880E-02 , 

4 1. 77085690E-02/ 

DATA ETA (01,14) , ETA ( 02 , 1 4 ) ,E TA ( 03 , 1 4 ) , E T A { 04 , 1 4 ) , ET A ( 05 ,1 4 } , 

* ETA( 06,14) ,ETA(07,14) , ETA ( 0 8 , 1 4) , E T A ( 0 9 , 1 4 ) ,ETA( 10,14) , 

* ETAl 11 ,14) ,ETA(12 ,14) ,ETA( 13,14) iETA{ 14,14) ,ETA{ 15,14)/ 

* 3.44827580E-02, 6 . 8 9 65 5 17 CE- 02 , 

1 9. 9616858 OE -02, ■ 1.226 05 360E-01 ,. 1 . 34865 900E-01 , 1 .3486 5900E-0 1 , 

2 1.23138430E-01, 1 .0 23488 20 E-0 1 , 7. 67616190E-02, 5. 1 1744120E-02 , 

3 2. 9627291 OE- 02, 1 .43647470E-02 , 5 .49240340 E-03 , 1 .47872400E-03, 

4 6 .8109622 OE— 03/ 

DATA ET A( 01 ,15 ) ,ETA (02,15) ,ETA( 03,15 ) ,ETA( 04, 15) , ETA C 05,15) , 

* ETAl 06,15) ,ETA{ 07,15) ,ETA( 08,15 ) ,ET A (09,15 ) , ETA (10,15) , 







o n 


* ETA(U,15),ETA(12,15) ,ETA{ 13, 15 ), ETA l 14, 15 >, ETA { 15, 15) / 

* 3.22500640E-02, 6.45161290E-02, 

1 9.3437 153 0E-02, 1 . 15 684090E-0 1, 1 . 2853 7S 80E- 01 , 1. 305 1 5390E -Ol , 

2 I. 2181 436 0E -01 , 1 - 044123 10 E~0 I s 8 . 17 139820E-02 , 5. 77775630E-02 , 

3 3.63L73250E-02, 1.98 0945 00E-02 , 9.03588950E- 03 , 3 . 2436 5270E-03 , 

4 8.17727570E-04/ . 

C ETA { I t J] 1=1,2,*.. , J IS USED IN THE FIRST MODIFICATION OF THE 

C I-TH DIFFERENCE OF A J-TH ORDER -METHOD AFTER THE STEPSIZE IS 

C HALVED. 

C ETA(I,J) J = 1,2,.*., .1-1 IS USED IN THE SECOND MODIFICATION OF 

C THE ( J+l )-S T DIFFERENCE OF AN I-TH ORDER METHOD*'-- - - 

C THE TWO MODIFICATIONS OF THE DIFFERENCE TABLE AFTER HALVING THE 

C ■ ' STEPSIZE REMOVES' MOST OF THE INSTABILITY INHERENT . IN THE METHOD.. 
C-- -USED HERE F OR-. HALVING -THE STEP S I ZE 

c . . : 

c - • - . 

C IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE FOLLOWING CARD. 

DATA LOSS ,LGSD,LGS E/0, 0,0/ . . .. 

INITIALIZE 

KSTEP=— 1 .... 

NE=NEQ 

IF (NE.LE.O) GO TO 1190 
HH=H 

NV = 0 - ' • 

KDMAX=0 . 

KDD=KD(1) 

KDS=KDD 
DO 10 1 = 1, NE 
KQ ( I ) =1 
DT{ 1 , I J =0* D 0 

IF (KDSdLEcO) kdd= I ABS { K D{ I ) ) 

IF { IKDD.EQ.OJ.OR. { KDD.GT.KMAXO) ) HH=O.DO . ......... 

IF ( KOD.GT .KDMAX) . K-OM AX=KDD - - . .. 

10 NV=NV+KDD ‘ ...... 

C ••••'. ‘ 

IF ( (DELT*HH) .LE.O-DO) GOTO 1190 
ERRMX=P1 . • '*• 

EMAX=ERND 
RNDC=RND*P2 5 

LDQUB=0 • • • 

E2HFAC = P 25 ■ 

LSC = 8 

LSTC=4 - 

C LSC AND LSTC ARE, USED IN COMBINATION AS FOLLOWS 


C 

LSTC =4 , 

LSC=4 

FIRST TIME THROUGH THE FIRST STEP 


c 

LSTC=3, 

LSC = 4 

SECOND TIME THROUGH THE FIRST STEP 


c 



(NECESSARY TO CHECK STABILITY) 


c 

L STC = 2 , 

L SC = 4 

THIRD TIME THROUGH THE. FIRST STEP 


_ c 



(ONLY OCCURS IF INSTABILITY POSSIBLE) 

c 

LSTC = 2 , 

LSC=2 

SECOND STEP (IF KQ(I)=2 , 1=1,..., 

NEQ ) 

c 

LSTC=1, 

LSC=0 

STARTING, ONE DERIVATIVE EVAL. PER 

STEP. 

-C 

. LS TC = 1 , 

LSC.GT.O 

SET WHEN STARTING TWO D E R IV. EVAL. 

PER STEP 

c 

LSTC =- 1 

LSC. LI. 0 

SET WHEN HALVING THE STEPSIZE 


c 

IN THE LAST TWO CASES LSC IS SET EQUAL TO L STC* l MAX IMUM 

KQ ( I ) 


C +1). AT THE END OF EACH STEP If LSC.NE.O IT IS REPLACED BY 

C LSC-LSTC UNTIL LSC=0, AT WHICH TIME LSTC IS SET EQUAL TO 0. 

C WHEN DOUBLING H, LSTC IS SET EQUAL TO -1 AND LSC TO -3. 

C UNDER CERTAIN CONDITIONS WHEN KQ(I)=1, LSTC IS SET =-l AND LSC=-5 

C 





KSOUT=MXSTEP 
TGUT=T 
I F L= 1 3 


20 I F LAG= 1 
GO TO 315 

C END OF I NIT IALIZATI ON 

C 

C ’ 

ENTRY DVDQ1 

- C 

c T ' _ 

c TO OUTPUT VARIABLES IN THE CALLING SEQUENCE REMOVE THE C 

* C -IN COLUMN ONE. OF -THE -FOLLOWING CARDS UNTIL REACHING THE 
C END OF CODE FOR PRINTING VARIABLES IN CALLING SEQUENCE. 

C IF CNEQ.NE.-O) GO TO 28 

C NEQ=1 

C 22 WRI T E! 6 , 500 0 ) T ,DE L T, HMI N A ,HM AX A , KE MAX , EMAX ? I FL AG « T F INAL 

C 5000 FORMAT l 3H0T= 1PD24. 1 7 , 7H DEL T=D 12. 5 , 8H HMI N A= ? E 1 0 . 3 , 8H 
C i El 0.3 * 8H KEMAX=*I2»7H EM AX= E10 . 3 , 8H IFIAG=,I2/ 

.C 2 9H I KQ KD,7X,4HF II > ,9X,1HJ ,12X ,4HY ( J) ? 22X,5HYN( J> , 

C 3 10X,7HTFINAL=1PD15.8 ,9H MXSTEP=I4} 

- C - - J=1 : 

C DO 24 X=1tNE 

C IF (KDS.LT.O) KDD-I ABS (KD( I) ) 

£ K = KDD 

C WRIT E { 6 ,500 1 ) I , KQ { I } , KDD , F { I)»JrY(J)»YN(J) 

C 5001 FORMAT ( 1 H , 12 ,21 3 , 1 PD1 7. 8 , 14 ,2D26 . 17 ) 


C 23 J=J+1 

£ K = K-*1 ... 

C IF (K.EQ.O) GO TO 24 

'C WRIT E ( 6 , 5002 } J»YU),YNU> 

C 5002 FORMAT! 26X , 14 » 1P2D26. 1 7 ) 

C GO TO 23 

*C ‘ 24 CONTINUE 
C WRITE! 6,5003) • 

C5003 FORMAT ( 3H0 I , 15X , 16 HD l FFER ENCE TABLE} 

C DO 2 7 I = .lt N£ 

. C KQQ=KQ ( I)+l 

C K=MI NO !KQQ» 7 ) 

C WRITE! 6,5004) I , (D T { I 0 , 1 ) 1 10=1 , K) 

-C5004 FORMAT ! 1H , 12 , 1 PD19 .8 , 6D 16 .7 ) 

C • ■ - ■ ■.-■■■ . ■■ ■ . • ’ .. - ■ .•• 

C IF ( K.EQ.KQQ) GO TO 27 

£ K = K+1 

C WRITE! 6,5005} (DT( 10,1 ),IO=K,KQQ) 

C5005 FORMAT! 1H , IP D2 1. 5 , 7D1 4 . 5 ) 

C 27 CONTINUE .... 

C IF (NEQ.EQ.O) RETURN 

C NEQ=0 

- C 28 CONTINUE 

C END OF CODE FOR PRINTING VARIABLES IN CALLING SEQUENCE. 

C 

IF I2-IFL) 30,60,320 
30 IF (IFL.GT.5) GO TO 1180 

C - 

. C SET STEP STOP 

. KSOUT=KSTEP +MXSTEP 
IF l I FL“4 ) 40,1210,210 
C 

C - 






-s 

COMMENT 


, MX STEP 
HM AX A= 
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C 


S E T PRINT STOP 
40 TOUT -T+DELT 


c : 

50 TPSI=ABS(SNGL (DMOD ( (TOUT-T )/ HH, 2 .DO ) J-PTSX) 

LF0=-1 

IF (TPS1.GE-P5) L FD= 1 

LFD IS USED TO INDICATE WHETHER DOUBLING H IS PERMITTED. 
IF LFD. L T. 0 AT THE END OF A STEP THEN DOUBLING H IS 
NOT PERMITTED. THE SIGN OF LFD IS CHANGED JUST BEFORE THE 
END OF EACH STEP. IF DELT#H*£POWER OF 2< THEN 
OUTPUT VALUES WILL BE OBTAINED WITHOUT INTERPOLATION. 

GO’ TO 200 1 


ENTRY WITH I FLAG=2 

* . • 

UPDATE DIFFERENCE TABLE - 

AND COMPUTE KQM-MAXI MUM VALUE OF KQU), I -l ,2 , . . . , NE Q. 

60 KQM=0 
• ' DO 80 1= L, N E 
KQQ=KQ( I ) 

IF ( KQQ. GT*KQM) KQM=KQQ 
D( l ) = F(I> 

DO 70 K = 1 » KQQ 
D( K+1)=DIK)-DT1K»I). 

70 DT { K * I > =D { K ) 

DT ( KQQ+1 ^ I ) -D IKQQ+1-) 

80- CONTINUE 

END OF UPDATING DIFFERENCE TABLE 

STORE Y( J) IN YNI J) 

DO 90 J-l t N V • 

90 YN { J ) = Y IJ ) 

LFD=~LFD ' ' ‘ • ' . > • ' 

TL-T • 

KSTE P=K S TE P +1 • • 

IF THF G STOP FEATURE IS ELIMINATED, REMOVE THE 2 FOLLOWING CARDS 
IF { L GS S ) 1430,110,1510 

100 I FL AG=2 ■ 

110 IF (LSC.EQ.OJ GO TO 140 . ' • 

LSC^LSC-LST C 
IF (LSC.EQ.O) GO TO 130 
IF ( LSTC *NE . { -1 ) J GO TO 140 
IF (LDOUB.LT. 0) RNDC=RND*P1 - 
120 E2HA VE-E 2HMAX 
TPSl^PTS 1 
GO TO 190 

130 IF ( A B S ( SN G L ( HH I) « L T • H M I N A ) GO TO 1000 

LST C=0 : 

140 IF ( LDOUB.NE* 1) GO TO 150 

IF { ( LFD.GT .0) .AND. CABS ( SNGL (HH+HH) ) •LE.HMAXA) ) GO TO 1030 
GO TO 200 

150 RQMA X=PT $1 /F LOAT l KQM+3) 

IF (LSTC.NE.O) GO TO 120 
TPS1=E2HMAX/£2HAVE ' • 
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IF (TPS1-PTSI) 160*190,170 
160 62HF AC = AMAX1 ( , 075EO , E2HF AC-RQM AX , E2HFAC*T PS 1 ) 

GO TO 180 
170 TPS1=TPS1*TPS1 

E2HFAC=AMIN1(PTS1»E2HFAC*TPS1) 

100 RN0C=( 1.1-E2HFAC) *RND 

E2HAVE=P5*i E2HMAX+E2HAVE) 

190 ERRMX=AMAX1( PI ,ERRMX-RQMAX*TPS1 } 

E2HFAC IS A FACTOR WHICH I S 'TAKEN TIMES AN INITIAL ESTIMATE OF 
E2H TO GET A FINAL VALUE OF E2H . IE2H=EST IMAT E OF WHAT 
{ESTIMATED ERROR I /( REQUESTED ERROR) WOULD BE IF H WERE 
DOUBLED* ), 

E2HMAX IS THE MAXIMUM VALUE OF THE INITIAL ESTIMATE OF E2H OVER., 

... ALL COMPONENTS WITH K0lI).GT.l. 

E2HAVE IS A WEIGHTED AVERAGE OF PAST VALUES OF E2HMAX. 

THE VALUE OF E2HFAC TENDS TO BE SMALLER WHEN E 2HMAX IS 
CONSI STANTL V SMALLER THAN E2HAV E . 


CHECK FOR PRINT STOP AND FOR T REACHING T FINAL. 

200 TPD=(TQUT-TL J/HH 

- TPD1=(TFINAL-TL)/HH 

IF THE GSTOP FEATURE IS ELIMINATED* REMOVE THE FOLLOWING CARD. 
IF (LGSE.LT.O) GO TO 1780 
IF (TPDl.LT.FACim GO TO 1220 
IF { TPD.LE. G. DO ) G'O TO 1280 

CHECK FOR STEP STOP 
IF { KS OUT .GT.KSTEP) GO TO 210 

I F L = 5 . • 

GO TO 310 ... . .. 

CHECK TO SEE IF ROUND-OFF ERROR IS PROMINENT 
210 IF ( EMAX *EQ • ERND) GO TO 220 

IT IS ....... - ' . \ 

I FL = 6 

IF (EMAX.GE.Pl 1 GO TO 310 

IF ( (LSTC.GE.O) .OR. (LDOUB.EQ.l) ) ERRNX=PTS1 

220 I FL = 1 • ; 

230 T=TL+HH 

START A NEW STEP ' 

PREDICT : 

. 2A0 J = 0 : . 

- DO 290 1=1, NE .. 

IF (KDS.LE.O) KDD=I ABS ( K D { I ) ) 

KDC-KDD 
250 KQQ=KO( I ) 

TPD=O.DO 

K=KDC 

260 TPD=TPD+DT( KQQ,I)*GAM{ KQQ , KDC ) 

KQQ=KQ0-1 

IF (KQQ.GT.O) GO TO 260 
. 270 K=K- 1 

IF (K.LE.O) GO .TO 280 • - 

L = J+ K 
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T PD-YN { L + 1 ) * F AC ( K )+HH*TPD 
GO TO 270 
280 J=J+1 

Y{ J)=YN( J ) + HH*TPD 

• - • KDC=KDC— 1 - 

IF ( KDC *GT . 0 ) GO TO 250 
290 CONTINUE 

END OF PREDICT . 

IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE C IN COLUMN ONE 
OF THE 2 FOLLOWING CARDS 
IF { IFL ) 20,320,300 
300 CONTINUE 

AND THEN REMOVE THE 2 FOLL OW I NG CARDS* :• , -■ 

IF { IFL ) , 1 240 ,320., 3-00 "" - ' . . . . . ? . " 

300 IF [ LGSD » NE * 0) ,G0 TO 1520 

310 IF LAG=I FL f • 

315 CONTINUE 

TO OUTPUT VARIABLES IN THE CALLING SEQUENCE REMOVE THE C IN 
; COLUMN ONE OF THE FOLLOWING CARD* 

IF ( NEQ. EQ* 0 ) GO TO 22 

. RETURN • ... 


ENTRY WITH IFLAG=1 
320 EPS=EP ( 1 ) 

ERND=0* 

EMAX=0* , . . ~ • 

E2HMAX=0. ' .... 

j=o . • . - ; ; •' ' •' 

IF (LDOUB.GE.O) LD0UB=1 ; - . 

LD0U8 IS SET IN THE LOOP BELOW AS FOLLOWS 

LDQU8=0 HALVE ...... 

LDOUB^l DOUBLE \ - 

LDOU B=2 DO NOT DOUBLE 

LDOUBy L T *0 AT THE BEGINNING OF THE LOOP I ND I CA TES THE FOLLOWING 
=-3 STEPSIZE HAS JUST BEEN HALVED. IF A DISCONTINUITY IS 
NOT INDICATED MODIFY THE DIFFERENCE TABLE AND REPEAT 
THE STEP. 

STEP AFTER LDQU0--3. PROCEED AS USUAL (ORDER IS NOT 
CHANGED) 

=-l STEP AFTER LD0UB=-2. MODIFY THE DIFFERENCE TABLE ONCE 
AGAIN AND REPEAT THE STEP. 

IF LDOUB IS SET EQUAL TO -4 THE ORDER IN AT LEAST ONE COMPONENT 
HAS BEEN GREATLY REDUCED AND THE STEP IS REPEATED. 


IF THE OUTPUT OPTION IS ELIMINATED, REMOVE THE 4:F0LLCWING .CAftDS 

IF . (NEQ. LE. 0) WRITE .(6,5020)., LSC,LFD,LSTC, KSTEP»E2HF AC ,ERRMX»HH 

5020 FORMAT (L9H0 I KQQ LRND L DO UB , 5X , 1HE ,9X ,3HE2H , 

1 8X,3HEPS,3X,4HLSC”» 13, 6H LFD=,r2,7H LSTC=,I2,8H KSTEP=,I4, 
•2 9H E2HFAC=,F4.2,8H ERRMX-»F4.2,4H H~,1PD9.2) 


BEGINNING OF LOOP FOR CORRECTING, ESTIMATING THE ERROR, 


: 'rry 5r , h 
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AND ADJUSTING THE NUMBER OF DIFFERENCES USED 
DO 790 1=1, NE 

IF { KDS , L E« 0 ) KDD= I ABS C K D ( I ) ) 

KQQ=KQ{ I ) 

KOQ GIVES THE ORDER OF THE PREDICTOR FORMULA AND KGQ+1 THE 
ORDER OF THE CORRECTOR FORMULA, 

KQ1=KQQ+ I 
D{ 1 ) =F ( I } 

FORM THE DIFFERENCE. TABLE FROM PREDICTED DERIVATIVE VALUES, 
DO 330 K=1 , KQ1 
D(Kf l)=D(K}-DT(K f I) • . 

330 -CONTINUE—: . - •• '. ‘ - ■ 

D(K) GIVES THE (K-l)-ST DIFFERENCE FORMED FROM PREDICTED 
DERIVATIVE VALUES 
TPS3=ABS (SNGLIDIKQQ+1)} ) 

I F { L00U3, L T- 0) GO TO 720 

340 IF (KQQ,NE. II GO TO 520 

KQ ( I )=1 IS TREATED AS A SPECIAL CASE 
E2H = P TS2 
TPS5=DT {3,1 ) 

IF (LSTC. LT ,2) GO TO 370. 

C FIRST STEP OF INTEGRATION 

IF ( LST C,NE ,4 ) GO TO 350 
TP S4=0, 

IF C KDD, GT , i ) TPS3=AMAXl(TPS3,ABS(SNGLlHH*OII)n ) * 

- ■ . TRS3-TPS3*P1 . . .• . 

. GO TO. 510 ‘ • 

. 350 DTI 2, I ) =01 2) , - . ■ • 

D(2)=D(U-DT(5,I) 

TP S2=~D { 2) - - 

TPS3=PTS5*ABSITPS2 J 
C FIRST STEP THAT KQ(I)=I 

- 360 DTI 7,1 ) =PT ( 4) 

IF (LSTC-2) 420,380,300 
370 IF ( TPS5 , EQ . 0 • ) GO TO 360 

IF (DT(6,I) ,EQoO,DO) GO TO 400 
TPS2=DT (5,1) — DT{I,I ) 

380 TP $4=DT( 4 , 1 ) 

TPS I=ABSITPS4) 

TPS4=TPS2*SIGN{PTS2»TPS4)— TP$5*TPSi 
IF ( TPS4.GT. (-TPSIU GO TO 410 
-- 390 TPS6=~PTS1 
GO TO 450 

C FIRST STEP AFTER THE STEPSIZE HAS BEEN CHANGED 

400 DT ( 6 , I )=PT ( 1) ... 

TPS 6=0. ■ : 

GO TO 450 ... 

410 IF (TPS4.LT.TPS1) GO TO 440 
IF ITPS1.E0,0.) GO TO 390 
420 TPS6=PTS1 
GO TO 450 
430 KQI I ) =2 

IF I2-LSTC) 510,510,520 
440 TPS6=TPS4/TPS1 

450 TPS4=TPS5*-TPS6 • 

IF (TPS4.LT ,P25} GO TO 430 
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C INCREASE E2H IF (-SKGT..25 

E2H=PTS4$TP S4 
IF { 2— L ST C ) 460,470,430 
460 L SC = G - " 

GO TO 510 

470 IF (TPS5-P25) 430,460,460 
480 IF (TPS4.GT *PTS2) GO TO 490 

IF (TPS4.GT.P5) .D(2) = D12)*GAM(2,1)- 
GO TO 510 

. . 490 IF { TP S4-LT*PTS4) GO TO 500 - ■ 

TPS4=PTS 4 
D(2)=D(2)/PT(3) . 

- C THE ESTIMATE OF E (AND HENCE OF E2H ) IS INCREASED IF {-S)*G£.8. 

. .. .. : -TPS3-TPS3*SNGL(DT.(7rI-) ) - : 

GO TO 510 

500 D{2)=D{2)*DBLE(PTS2*(TPS4-PTS1)/(TPS4*TPS4} ) 

IF ( TP S4 oGE *P3E 1 ) E 2H=E 2H* SNGL ( DT { 7 , I )) 

STORE Dll) ^PREDICTED DERIVATIVE AND DC 2)- 2* (CORRECTED Y- 
PREDICTED Y)/H Dll) AND D(2) ARE USED TO COMPUTE (-S) ON 
THE NEXT STEP. 

510 DT ( 5 , I ) =D( 1 ) 

. DT ( 4, 1 )-Di 2) 

D( 4 )~TPS4 

' STORE D 1 4) = CURRENT ESTIMATE OF (-SK (-SKGT.3 IS AN INDICATION 
THAT THE STEPSIZE SHOULD BE LIMITED BECAUSE OF STABILITY PROBLEMS. 
S=H* { EST I MAT £ OF EIGENVALUE OF F) = H* ( D IFFERENCE BETWEEN PREDICTED 
AND CORRECTED DERIVATIVE VALUE S )/( DIFFERENCE BETWEEN PREDICTED 
AND CORRECTED INTEGRALS OF THE DERIVATIVE VALUES) 

THE TREATMENT OF THE CASE KQII)=1 COULD BE IMPROVED BY USING A 
SPECIAL METHOD FOR STIFF EQUATIONS WHEN (-S)„GT.3 (MAYBE). 

(THE ENTIRE TREATMENT OF THE CASE KQ(I)=1 IS FAR FROM IDEAL.) 

DT (3,1 ) = D ( 4 ) 

CORRECT ’ ’ . ..... 

' 520 KDC = 0 ..... .. .... . .. 

T PD=D ( K Q1 ) ■ - ; 

J= J+KDD ... • ••.- 

K=J 

53 0 T l PD=HH*TPD ... . 

KDC=KDC+1 

. Y(K)-Y(K)+GAM(KQQ+1 ,KDC)*TPD 

K= K~1 

IF (KDC.LT.KDD) GO TO 530 

END OF CORRECT . . 

IF ( EPS ) 540,550,560 
540 EP S = EP { I ) 

IF ( EPS *NE • 0 • ) GO TO 560 
550 IF IHMAXA) 1190,780,1190 
560 TPS4=ABS(SNGL(D(KQQ+2) ) ) 

TP52=ABS (SNGL (D(KQQ )) ) 

TPS6=SNGL(HH) /EPS 

E=ABS (SNGL ( GAS(KQQ+ 1) )*TP$3*TPS6) 

E GIVES AB S( (ESTIMATED ERROR)/EPS) 

LR ND-1 ’ 

LRND= 1 MEANS NO ROUND-OFF ERROR 
« 0 MEANS SOME ROUND-OFF ERROR 



c 

•c 


=-i means extreme round-off error 


B20 


FRND=RNDC*ABS(SNGL ( PT [ KQQ+2 ) *D{ 1 ) ) ) 

C CHECK TO SEE IF ROUND OFF ERROR IS DOMINANT 

IF ( (TPS3+TPS4.)— GT-FRND) GO TO 570. 

LRND=0 

IF ( (PTS4*TPS2ULT.FRND> LRMD=- 1 
C ■ ■ 

570 IF (E a LE.ERNDJ GO TO 590 
IF fE.LE.EMAX) GO TO 580 

emax=e 

KEMA X= I - 

580 IF (LRND.LE.OJ GO TO 590 
ERND=E 

IF (ERNDeGT.ERRMX) LDOUB^O 
• 590 IF (LOOUB* LE.O) GO TO 780 

TPS l =A SSI SNGL IDD(KQQ)) ) * 

TPS 5 -TPS 1 

IF'IKQQ-2) 600,610, 620 
-- 600 E2H=E*E2H ... 

IF I E2H * LT a P01 ) GO TO 780 

IF ( SNGL ID ( 4J ) .LT . P3EI ) GO TO 770 • * 

LSTC=-1 . ■ 

LSC=-5 . ... ... 

>-:■ - GO. TO 770 ■ . .. ... 

610 TPS1=TPS2 ‘ * 

IF (LSTC.NE.2) GO TO 620 
KQ ( I ) =3 

TP S2 = 0. • 

TPS 4-0 « 

. LRND=0 

620 £2H=TPS2 + TP 53+TPS4 

E2B=ABS (SNGL! GAS (KQQ-1 } *PT { K QQ+ 1 ) )*E2H*TPS6) 

'C E2H IS USED AS AN ESTIMATE OF WHAT THE VALUE OF E WOULD BE 

C IF H WERE DOUBLED. THE ESTIMATE IS CONSERVATIVELY LARGE. 

IF ( E2H.GT * E2HMAX ) E2HMAX-E2H 

G ■ • . ’ 

' IF { LRND ) 630,640,660 
C EXTREME ROUND-OFF ERROR — REDUCE E2H 

630 K=IKBIT2/KQQ) —4 

IF (K.LE.3) GO TO. 640 
, IF . (K.GT.KQMAX) K^KQMAX 

E2 H= E2H/PT { K+l) ... 

GO TO 650 

640 E2H=AMrNllE2H,E2H*P3EI*e2HFAC) . - 

, 650 E2H=E2H*P1 

TPS6-PTS4 
GO TO 670 

. C ... 

660 E2H=E2H*E2HFAC 

TPS 6= FLOAT { KQQ+-2) 

c TEST TO SEE IF .DIFFERENCES DECREASE MORE RAPIDLY THAN. NECESSARY 

-C 

670 IF ( TPS 5.LT . { TPS3*TPS6) } GO TO 680 
IF (TPS2.LE. (TPS4*TPS6 ) } GO TO 760 
C THEY DO INCREASE KQ(I) 

IF ( KQQ.NE.KOMAX) KQ(I)=KQI 
GO TO 760 
C 

C TEST TO SEE IF DIFFERENCES DECREASE TOG SLOWLY 
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680 TPS6=TPS6*.P25 

IF ( (TPSi.GT. (TPS3*TP$6) ) * OR . ( T PS2 . G T* (TPS4*TPS6 ) ) ) GO TO 760 
C THEY 00 

IF (LSTC.LE.O) GO TO 750 . . 

IF { E2H *LT * PO 1 ) GO TO 750 
IF { LSC-LST C) 690,750,770 
. 690 IF (KSTEP-4) 750,700,710 

700 KQ1 =LST C - . . 

710 LSC=KQ1 . .. 

- C END OF ONE DERIVATIVE EVALUATION PER STEP 

- • GO TO 770 ; 

AFTER HALVING H* REDUCE KQ ( I > IF A DISCONTINUITY HAS OCCURRED* 
720 IF ( LDOUB *EQ- (-2) ) GO TO 340 

DT (KQQ+1 »I)=D(KQQ+1 ) . 

IF (LDOUB.EQ. (-1) > DT ( KQQ+1 , I ) =D { K QQ*2 ) 

K=KQQ 

730 IFiK.EQ.l) GO TO 740 ... 

IF ( (DABS(D (K-l) ) eGT.I PT( 2)*DABSID{ K+l) M } . OR* 

1 ( DABS ( D { K ) ) *GT * ( PT ( 2 ) *DA BS(D(K+2))))) GO TO 740 

K-K- 1 

GO TO 730 . .. . 

740 IF ( (K+Kl.GE.KQQ) GO TO 780 . 

LD0UB=~4 . 

E2H=0 • 

KQQ=K+1 


DIFFERENCES DECREASE TOO SLOWLY REDUCE KOI l I * 

■750 KGU)=KQQ-1 . ...... •• 

• IF • ( KQQ . EQ • 2 ) DTC 3 ,11 -0*D0 .. .. 

760 IF [ E2H.LT.P01 ) GO TO 780 
770 LD0UB=2 

780 CONTINUE . 

IF- THE OUTPUT OPTION IS ELIMINATED, REMOVE THE 6 FOLLOWING CARDS * 
IF { N £ Q . GT .0 ) GO TO 790 
TO2 = MAX0U, CKCQ-1) ) 

103=102+3 

WRITE {6 ,50211 I , KQ Q, L RND , LDOU B , E , E2H , EPS , 

1 UG1,D(I01) ,101=102,103) 

5021 FORMAT ( 1H 3 2 , I 4, 2 1 5 , IP E 13 . 3 , 2E 1 1 . 3 , 4{ 3H { , 1 2 , 1H ) , D 1 0. 3 )) 

C 

790 CONTINUE 

C - 

.C END OF LOOP FCR CORRECTING, ESTIMATING THE ERROR, ETC* 

C 

C 

C IF THE INTERPOLATION CAPABILITY IS ELIMINATED REMOVE THE 

-C FOLLOWING CARD. , 

I.F ( XFL *LT*0). GO. TO .1250 . - . 

C-«* TEST FOR- HALV ING--H , . 

" IF ( LDOUB) 80 C, 95 0, 870 

800 LDOUB=LDOUB+i 

IF (LDOUB+1) 810, 870,820 - - : 

810 IF (LDOUB. EQ. (-2) ) GO TO 820 
C ORDER IN AT LEAST ONE COMPONENT HAS BEEN GREATLY REDUCED 

LD0UB=0 
GO TO 220 
820 DO 860 1=1 ,NE 


wry? 7 ’" 
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KQQ-KQ( I ) 

TP = DT ( KQO+1 , I ) 

IF (KQCULE.3) GO TO 860 
IF ( LDOUB * NE . 0) CO TO 840 
DO 830 K=3 » KQQ 

C SECOND MODIFICATION OF DIFFERENCE TABLE AFTER HALVING H 

83 0 DT (K, I ) = DT ( K, I) + ET A (KQQ- 1» K-2) *TP 

GO TO 860 .■■■•-- ' ■ 

840 OD 850 K = 2 , KQ Q 

- C FIRST MODIFICATION OF DIFFERENCE TABLE AFTER HALVING H 

850 DTlK,l)=DT(K,I)+ETA(K-l,KQQ-l)*TP 

■ , 860 CONT INUE- • . . ....... . ... . . \ 

IFL=0 ■ 

GO TO 240 
C 

870 I F L~2 ... 

IF (LSTC.LE *0} GO TO 300 
IF (2-LSTC) 880,900,940 . . 

• • 880 L STC = L S TC-1 

IF (LSTC.EQ.3) GO TO 890 
IF (LSC) 920*960,920 
890 IFL=I 

GO TO 300 - ■ 

•900 IF ( LSC-2) 910,930,920 ‘ 

• • 910 L S TC = 0 

920 L DOU B=2 

GO TO 60 

- ■ 930 LS TC = 1 

LSC=0 

GO TO 60 

940 IF (LSC) 300,60,300 
HALVE H 

950 HH=FAC(2 )*HH 

IF (LSTC.LT-2) GO TO 990 
ERND=P25*ERND 

IN LOOP TO FIND A NEW INITIAL STEPSIZE 
IF (ERND.GE.Pl) GO TO 950 

L S T C = 4 

960 LSC = 4 

DO 970 1=1, NE 

970 KQ ( I )"1 .... 

IF (LSTC-3) 890,890,1170 

ENTRY AFTER I FL AG=7 
9 BO IF (LDOUB. EQ.O) GO TO 990 
LS C= 1 
LSTC=l 
GO TO 140 

- C TEST TO SEE IF H IS TOO SMALL FOR HALVING 

990 IF ( A B S ( SNGL(HH) ) .GE.HMINA) GO TO 1040 
IF ( IFL.EQ.7) GO TO 1010 
1000 TFL=7 

GO TO 1020 

C - 

1010 HH=HH+HH 
l F L = 2 
1020 H=HH 

GO TO 310 .... 


C 
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ERROR CRITERIA PERMIT DOUBLING 
1030 HH=HH+HH 

IF (LSTCoEQ-1 ) GO TO 1050 
LSC=-3 

1040 LS T C=- 1 

CHANGE THE STEPSIZE ’ ,, 

1050 DO 1160 1=1, NE 
KQQ=KQ { I ) 

IF ( KOQoNEo 1) GO TO 1070 ' • 

DT { 6 s I ) =0 o DO 
D 1 3 ) =DT ( 3 » I )*PT{2) 

IF- ID( 3) .GT.PT<3)) LSC=-6 
IF . (LDOUB«NE*0) GO TO 1060 
KQM-8 

IF ( D 1 3 ) -GE *PT {5 ) ) DT ( 7 , I )=DT{ 7 , I > *PTl 2) 

D(3)^D(3)/PT(3) - 

1050 DT(3fI)=D{3) 

GO TO 1160 

BEGINNING OF LOOP FOR CHANGING DIFFERENCE TABLE TO 
CORRESPOND TO NEW VALUE OF H 
1070 DO 1080 K=l, KCQ .... 

•D(K)=DT IK, I )/ PT (K) 

1080 IF ( LDOUB© E Q» 0) D { K } =0 { K ) / PT ( K ) 

KQQ2=KQQ-2 

IF IKQQ2) 1 160 , 1 140 , 1090 

1090 DO 1130 J=1,KGQ2 ; 

IF (LDOUB.NE.O) GO TO 1110 

HALVE ; 

K = KQQ 

1100 D { K— 1 ) = D { K— 1 ) +D { K ) 

— K ■*" 1 - ■■.- 

IF (K+J-KQQ) 1130,1130,1100 

DOUBLE - " 

1110 DO 1120 K=J,KGQ2 
1120 D{K+l)=D(K+l)-D(K+2) 

1130 CONTINUE 

1140 DO 1150 K=2 , KQQ 

IF (LDOUB.NE.O) D l K )=D ( K)*PT ( K) 

DT(K,I)=D(K)*PT(K) 

1150 CONTINUE 

DIFFERENCE TABLE NOW CORRESPONDS TO NEW VALUE OF H 

1160 CONTINUE 
1170 H=HH 

IF ( LDGUB* NE* 0 ) GO TO 50 
LFD= 1 

IF ( LST C - GE *0 ) GO TO 220 
LOOUB =- 3 
LSC = LST C-KOM 
GO TO 220 

END OF CHANGING STEPSIZE 


IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE C IN COLUMN ONE 
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C OF THE 2 FOLLOWING CARDS 

*C1180 IF { 7- I F L } 1181,980,220 

C L 18 1 IF (IFL-8) 60,1200,60 

C AND THEN REMOVE THE 2 FOLLOWING CARDS. 

1180 K- 1 F L- 5 

GO TO (220,980,1200,1570,1570,1720,1720,60,1480,1450,1630,1570), K 
C 

C ILLEGAL VALUE OF PARAMETER INTEGRATION CAN NOT PROCEED 

1190 IFL-8 • - : '• 

• GO TO 310 

1200 VRITE ( 6,4000) 

4000 FORMAT (26H0IFLAG=8 IN CALL TO DVDQ1.) 

STOP . . 

C ' • 

C 

1210 IF (T-TFINAL) 200,1190,200 
C 

C . IF ONE DOES NOT WANT THE I NT ER POLAT I CN FEATURE, REMOVE ALL CARDS 
. C- • BELOW THIS POINT -{ EXCEPT FOR . THE END STATEMENT), AND ADD THE 
C FIVE FOLLOWING STATEMENTS* 

C 1220 I F L = 4 

C . IF (TPD1.GT .TPD) GQ TO 1280 

C GO TO 310 

C 1280 IFL=3 
. C GO TO 310 

• C 

1220 IFL=4 

IF { KS T EP. N E • 0 } GO TO 1270 
T P D2 = TPD • 

C ESTIMATE ERROR WHEN EXTRAPOLATION FROM INITIAL POINT IS REQUESTED 

.‘--1230 HH=HH*TPD1*.75D0 
C - - • 

..C . IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE FOLLOWING CARD. 

I FLS - 1 F L 

I FL = - 1 ... . ' . 

GO TO 230 

C . 

-C IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE 4 FOLLOWING CARDS. 

-1240 IF ( (LGSD.EQ.O) .CR.{ IFLS.NE.A) ) GO TO .20 . .. 

LG SE=- 1 

TP D= FAC ( 1) 

GO TO 1820 . * 

1250 HH=H 1 . 

IF (EMAX.LT.P01) GO TO 1260 ' . 

C ERROR IS TOO LARGE, REDUCE H. AND REPEAT THE FIRST STEP 

IF (TPDl.LT.O.DO) GO TO 1190 
LD0UB=1 

ER ND=F AC ( 1 ) /TPD1 
ERND=ERND*ERND*P25 . . 

- GO TO 950 

C 

C IF THE GSTOP FEATURE IS ELIMINATED, REMOVE THE C IN COLUMN ONE 

" C OF THE FOLLOWING CARD 

C 1 260 I F L=4 

C AND THEN REMOVE THE 2 FOLLOWING CARDS. 

1260 I F L= I F L S 

IF { IFL.NE.4) GO TO 1790 
TP D=TPD2 

• • - IFLAG= 3 . ; . 

1-270 IF (TPD1.GT .TPDl GO TO 1280 
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t=tfinal 

TPD=TP D1 
GO TO 1290 
12B0 T- TOUT 
I F L = 3 

1290 IF UTPO.EQ. 0.00). AND. { I F L AG „ LE . 2 ) ) GO TO 310 
C 

C INTERPOLATE FOR OUTPUT - ' ‘ ' 

1300 TP = TPD 
D{ 2 )=TP 

K0Q2 =0 ' •- 

KDC = 0 

Of 1 l = PT Cl) 

DD(i)=PT(l) 

DO 1310 K=2 ,KQM 
DD( 1 )=DD l 1 ) +PT( 1) 

TP = TP*PT (1) 

1310 D(K + 1)=(D(K )*TP)/DD{1) 

GO TO 13 50 V 

c • : 

C COMPUTE THE INTERPOLATING INTEGRATION COEFFICIENTS 

1320. KQG2 =1 

L= KQM- KDC ........ 

KDC=KDC + 1 . .. 

1330 IF (L.LE.O ) GO TO- 1350 * 

TP = 0. DO 
K = L • 

J= L+KDC 
1340 J S = J~K 

TP = TP+GAS(K )*D( JS-t-U *-• 

K=K~l 

IF (K.GT.O) GO TO 1340 . . 

D( J )-T P 

C - ' 

C D { J ) IS THE INTEGRATION CO EF F IC 1 ENT . FOR THE INTERPOLATION WHICH 

C CORRESPONDS TO GAM ( J-KDC , K OC ) . 

C - ....... 

L = L-1 

GO TO 1330 

C END OF COMPUTING INTEGRATION COEFFICIENTS 

C • • 

C PERFORM THE PARTIAL STEP INTEGRATION 

1350 J = 0 - - ■ 

DO 1420 1=1 ,NE 

IF (KDS.LE.O) KDD=IABS (KDC I) ) 

IF (KDC.GT.KDD) GO TO 1410 
TP = 0 .DO 

KQQ = KQ( I ) + K QQ2 
1360 L=KQO-KDC 

IF (L.LE.O) GO TO 1370 
TP=TP+D(KQQ)*DT(L, I ) 

KQQ=KQQ- 1 

IF { KQQ ) 1390,1390, 1360 
1370 K= J+KDD 
l=kdc 
1380 L=L-1 

IF ( L.EQ.O) GO TO 1400 * 

TP = TP*HH+YN IK )*FAC ( L)*TPD 

K=K-1 - : 

GO TO 1380 
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1390 F l I ) =T P 

GO TO 1420 

1400 Y{ K)=YN{K)+HH«TP 

1410 J=J+KDD " 

1420 CONTINUE 

IF ( KOC .NE • KDMAX) GO TO 1320 
END OF PARTIAL STEP INTEGRATION 

IF THE GST 0 P FEATURE IS ELIMINATED* REMOVE THE C IN COLUMN ONE 
OF THE FOLLOWING CARD 
GO TO 310 

ALL STATEMENTS BELOW THIS POINT SHOULD THEN BE REMOVED (EXCEPT 

~ FOR THE END STATEMENT) 

IF ( LGSE ) 1 800, 310? 181 0 


SECTION FOR COMPUTING GSTOPS ; 

ENTRY DVDOG l NG* NSTOP *G *GT ) 

C - 

C VARIABLES IN THE CALLING SEQUENCE- HAVE THE FOLLOWING TYPES. 

. INTEGER NG * NSTOP 
DOUBLE PRECISION Gtl)*GT(l) - 
C 

.. c A GSTOP IS DEFINED AS A RETURN WHICH IS MADE TO THE USER WHEN A 

C- USER SPECIFIED FUNCTION G PASSES THROUGH ZERO. THE USER MAY 
C SPECIFY ANY NUMBER OF FUNCTIONS G OF TWO TYPES * ZEROS OF THE FIRST 

C TYPE ARE LOCATED WITHOUT REQUIRING A DERIVATIVE EVALUATION 

C BEYOND THE ZERO. THIS TYPE OF GSTOP REQUIRES THAT G BE EVALUATED 

C .. BEFORE EACH DERIVATIVE EVALUATION. ZEROS OF THE SECOND TYPE ARE 

:.C LOCATED USING INTERPOLATION* WHICH IS MORE ACCURATE THAN THE 

C - EXTRAPOLATION USED IN THE PRECEDING CASE AND ONLY REQUIRES ONE 

C EVALUATION OF G PER STEP. THUS ONE SHOULD USE THE SECOND TYPE OF 

. *C GSTOP IF POSSIBLE. USERS NOT USING THE GSTOP FEATURE NEED READ 

C NO FURTHER. 

C ' 

C DVDQG IS USED AS A SET UP CALL TO INDICATE A CHANGE IN THE NUMBER 

. C - OR TYPES OF GSTOPS. DVDQG SHOULD BE CALLED JUST BEFORE OR JUST 

C AFTER CALLING DVDQ IF 

C 1. ONE WANTS TO TEST FOR GSTOPS BEGINNING WITH THE FIRST STEP. 

C 2. A JOB IS BEING RUN AFTER ANOTHER JOB THAT USES THE GSTOP 

C -FEATURE. DVDQG MUST BE CALLED EVEN IF ALL THE VARIABLES IN 

C THE NEW JOB ARE THE SAME. 

C IN ADDITION DVDQG MAY BE CALLED AT ANY TIME IN THE INTEGRATION 

C TO. CHANGE THE NUMBER OR TYPE OF GSTOPS. 

C 

C THE USAGE OF THE VARIABLES IS GIVEN BELOW. 

C 

C . NG = 

** C 

c 

C 

* G - 

C 

c 
c 
c 
c 
c 
c 


THE NUMBER OF COMPONENTS IN G TO BE EXAMINED FOR A ZERO. 

IF DVDQG IS CALLED AFTER THE FIRST STEP OF THE INTEGRATION* 
THEN G IS EVALUATED FOR THE FIRST TIME AT THE END OF THE 
NEXT STEP AND THUS A GSTOP IS NOT DETECTED IF G CHANGES. 

SIGN ON THE . CURRENT STEP.. IF .IT -IS .IMPORTANT THAT G BE 
EVALUATED IMMEDIATELY SETNG EQUAL TO THE NEGATIVE OF THE 
NUMBER CF COMPONENTS TO BE TESTED FOR A ZERO. SETTING NG 
LESS THAN ZERO WHEN CALLING DVDQG BEFORE THE FIRST STEP IS 
NOT NECESSARY AND IS LIABLE TO BE DISASTEROUS. IF DVDQG IS 
CALLED DURING THE INTEGRATION THE FOLLOWING STATEMENT SHOULD 
BE A GO TO (THE COMPUTED GO TO FOLLOWING THE CALL TO DVDQ1 ) . 




c 

■c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c- 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c- 
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NSTOP=THF NUMBER OF COMPONENTS OF G THAT MUST BE EXAMINED FOR 
A ZERO BEFORE COMPUTING THE DERIVATIVES (FIRST TYPE OF 
GS TOP ) - IF NSTOP.LT. 0 OR NSTOP . GT . ABS ( NG ) IFLAG IS SET 
EQUAL 8 AND AN IMMEDIATE RETURN IS MADE. IF NSTOP. GT.O, 

G(1 ) , G( 2 ) , . ..yG(NSTOP) ARE EXAMINED FOR A ZERO BEFORE EACH 
DERIVATIVE EVALUATION, THE REMAINING COMPONENTS (IF ANY) 

ARE EXAMINED AT THE END OF EACH STEP, WHEN A GSTOP IS FOUND 
THE SUBROUTINE SETS - NSTOP EQUAL TO THE INDEX OF THE 
COMPONENT OF G RESPONSIBLE FOR THE STOP. 

A VECTOR CONTAINING THE CURRENT VALUES OF THE FUNCTIONS 
WHOSE ZEROS ARE TO BE DETERMINED. 

GT = A VECTOR WITH THE SAME DIMENSION AS G USED BY THE 
SUBROUTINE FOR TEMPORARY STORAGE. 

RETURNS FROM CALLING DVDQ1 WITH I FLAG. GT *8 SHOULD BE INTERPETED 

AS FOLLOWS. (WE USE NSTOPI TO DENOTE THE INITIAL VALUE OF NSTOP.) 

IFLAG 

- 9 COMPUTE G(NST0PI+1) GIABS(NG) ) (THE COMPONENTS OF G WITH 
ZEROS TO BE LOCATED USING INTERPOLATION). THEN CALL DVDQ1 . 

NO RETURN IS MADE WITH I F L A G= 9 IF N STOP I = A B S ( NG ) . 

=10 COMPUTE G<1) ,G(2) ,... ,G(NSTQPI) (THE COMPONENTS OF G WITH 
. ZEROS TO BE LOCATED USING EXTRAPOLATION). THEN CALL DVDQ1 • 

NO RETURN IS MADE WITH IFLAG=10 IF NSTQPI=0. 

=11 G( NSTOP ) IS APPROXIMATELY ZERO. IF THERE ARE NO 

DISCONT 1NUITIES SIMPLY CALL DVDQ1 TO CONTINUE THE INTEGRATION* 

=12 G ( NSTOP ) CHANGES SIGN, BUT THERE IS DIFFICULTY IN CONVERGING 
TO A ZERO. THE -USER HAY WISH TO MAKE. A SPECIAL CHECK TO BE 
CERTAIN THAT EVERYTHING IS ALL RIGHT. TO CONTINUE THE 
INTEGRATION CALL DVDQ1 . . . . 

DOUBLE PRECISION RG 

DOUBLE PREC ISION GI . .. 

DIMENSION GI(2),RG(3) .. 


INITIALIZE FOR GSTOPS 
NGA=IA8S( MG ) 

LGSS=-NGA 
LG SD = 0 
LG SE = 0 
I F LG=— 2 0 


IF (NG) 1425,315,315 

• 1425 IFLG=-IFL 

I FLG=- I FL ' ■ ' 

1430 LG $D=N$T OP 

IF ( LGSD } 1 190,1450,1440 
1440 I FL=1 5 

GO TO 1470 

C ENTRY WITH IFL=15 

1450 LGS S=0 

IF ( LGSD-NG A ) 1460,1480,1190 
1460 LGS S = LG$D+1 
I F L= 14 

.. 1470 IFL AG= I FL- 5 

• •- GO TO 31 5 ... • .... 

C ENTRY WITH 1FL=14 

- 1480 DO 1490 1 = 1 ,NGA 


1490 GT ( I ) =G ( I ) 
GO TO 1730 
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C END OF INIT IALIZATIGN FOR GSTOPS 

*C 

C ENTRY TO EVALUATE. G AT THE END OF THE STEP 

1500 LGSE=1 
1510 I GK = LG S5 
I FLG=0 
I FL = 9 

GO TD 310 . . - 

C ENTRY TO EVALUATE G BEFORE EVALUATING THE DERIVATIVES 

- 1520 IFLG=IFL 

ifl=io 

1530 IFLAG=10 
IGkM=LGSD 
■ 1540 IGK=1 

1550 GO TO 315 ......... , 

.*-1560 I GK= I GK + 1 - ~ ... - , 

IF { IGK.GT. IGKM) GO TO 1650 

C ENTRY WITH IFL=9tlO» AND 17 . - * 

C TEST FOR G CHANGING SIGN 

■ • 1570 IF ( G{ IGK)*GT( IGK) ) 1600 » 1 58 0 , 1 590 

1580 IF ( GT { IGK) .NE* 0* DO ) GO TO 1600- 
IF (TL«EG.TG) GO TO 1560 

1590 IF (LGSE.GT.O) GT ( I GK ) =GU GK ) 

GO TO 1560 . 

C G CHANGES SIGN — PREPARE FOR ITERATION TO FIND ZERO 

-1600 NS TOP = I GK - -• 

NSTO PI = I GK ....... 

. .. I FLGS = I FL . , •' 

C COMPUTE INITIAL VALUE FOR RG (=RATIQ OF PARTIAL STEPSIZE WHERE 

.C G IS KNOWN/THE INTEGRATION STEPSIZE) 

. IF { I FLG * EQ .0 ) GO TO 1610 
RG( 3 ) -F AC { 1 ) 

* . RG(2)=0 -DO 

IF ( ( IFLG.E C.2) -AND. UGK.LT.LGSS ) ) RG(2) = FACU) 

GO TO 1620 , 

1610 RG ( 3 ) = 0- DO • • 

RG { 2 ) =- F AC l 1 ) •• 

1620 IF (LGSE-LT.O) RG(3)=TPD 

LGS E=-3 - 

GI (2 ) = GT( IGK) 

EP SGS=RND • ...... • : 

I FL=1 6 ; 

K= 1 . 

GO TO 1640 

C END OF PREPARATION TO BEGIN THE ITERATION . ■. 

C ENTRY WITH IFL=16 .. . ■ 

C ITERATE TO FIND GSTOP .* 

1630 K=1 

IF ({GI l 2 ) *G ( IGK) ).GT*0*D0) K=2 
IF l DABS ( GI ( K ) ) .GT * DABS ( G ( IGK ) ) ) GO TO 1640 
*_-C CONVERGENCE PROBLEMS * ■ 

LGSE=LGSE-1 

IF (LCSE.EQ.(-5H EPSGS=PTS1 
EPSGS=EPSGS*PTS4 
1640 GI ( K ) = G ( IGK ) 

RG ( K) =RG ( 3 ) 

C SECANT ITERATION (GIVES. NEW PARTIAL STEPSIZE/H) 

T PD = RG (1 )-( GI (1)*(RG12)-RG(1> ) ) / l G It 2 )-G I ( 1 ) ) 

T=TL+TPD*HH 
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TEST FOR CONVERGENCE OF ITERATION 
IF (DAB${TPD-RG(3) ) .LE * EPS GS ) GOTO 1560 
RG( 3 )=TPD 
GO TO 1300 

1650 IF (10-IFL) 1660,1700,100 

1660 IF { I GKMo NE * NGA ) GO TO 1710 
IF .(LGSE.GT.I-3J) GO TO 1690 
IF (LSTC.NE.4) GO TO 1670 
C ESTIMATE ERROR — G STOP IS T HE RESULT OF EXTRAPOLATING FROM 

C THE INITIAL POINT ' , 

TP D1 =T PD 
RG ( 3) =TPD 
GO TO 1230 
1670 I FL=1 1 

IF ( LGSE .LT . ( -4 ) ) I FL=I2 
1680 I FLAG=I FL 

C TEST TO SEE IF GSTOP IS PRECEDED BY . ANOTHER STOP 

IF { ( {T-TQUT)*HH.LE,O.DO) .AND. { ( T-TF I NAL ) *HH . LE . 0 . DO ) ) GO TO 1300 
C IT IS 

RG ( 3 ) =TPD 
I FLS = I FL 
GO TO 200 
1690 LGS E= l 

I FL= I FLG 

IF (IFL.LT.O) GO TO 20 
1700 I GKM=NGA 

• ■ IFL= TFLG - : . .. 

GO TO 310 
1710 I FL= 17 
I FLAG=9 

- - * I GKM=NGA " • •. • 

GO TO 315 : ... 

C ENTRY WITH I F L= 11 AND 12 

C SET PARAMETERS TO INOICATE A GSTOP HAS BEEN FOUND 

1720 GT (NSTOPI ) = 0 • DO ' 

1730 LG SE =1 - . 

IGKM~NGA 

TG=TL 

IF { TFLG ) 1740,1760,1770 
-1740 IF (IFL.LT.13) GO TO 1750 

IF (IFLG.EO. (-20) ) GO TO 100 
... . IFL=-IFLG ■ ■ 

GO TO 310 
1750 HH=H 

GO TO 200 , 

1760 T PD= 0 . DO 

-• - . T=TL - 

LG SE =- 2 
GO TO 1300 

1770 IF (IFLG-3) 220,200,200 
1780 IF (LGSE.EQ.(-l)) GO TO 1790 

LGSE^-l • ' •. 

* GO TO 1220 
1790 TPD=RG(3) 

. T=TL+TPD*HH . . - - - .. - 

IF (LGSE.NE .(-I)) GO TO 1670 

IFL=I FLS • . 

LG SE =-3 
G3 TO 1680 

1000 IF (LGSE+2) 1550,1500,310 
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1810 IF (TPD.LE.0.D0) GO TO 310 
LGSE--2 
1820 1FLG=IFL 
I F L = 1 7 
I FL AG=9 

IF ILGSD .GT. 0) GO TO .1530 
GO TO 1540 

C END OF SECTION FOR COMPUTING GSTOPS 

C 

c 

C IN SOME APPLICATIONS* FOR EXAMPLE MULTIPLE QUADRATURE* MORE THAN 

C ONE INTEGRATION SUBROUTINE IS REQUIRED. THIS IS NOT NECESSARY IF 

C ALL OF THE VARIABLES ASSOCIATED WITH ONE INTEGRATION ARE SAVED 

C OUTSIDE OF THE INTEGRATOR WHILE DOING OTHER INTEGRATIONS* AND.. . . 

G" * , THEN .RESTORING- THEM. WHEN. NEEDED.. THESE VARIABLES CAN BE RESTORED 
C BY CALLING -AN ENTRY WHICH CONTAINS ALL OF THE VARIABLES IN THE 

C CALLING SEQUENCE AND ALL OF THE VARIABLES THAT MUST BE SAVED 

C WHENEVER CHANGING TO A DIFFERENT INTEGRATION. ITHIS ENTRY MUST 

C BE ADDED BY THE USER AND SHOULD BE FOLLOWED BY A RETURN STATEMENT. 

C AFTER CALLING THIS ENTRY EITHER DVDQ OR DVUQ1 SHOULD BE CALLED 

C DEPENDING ON WHETHER THE INTEGRATION IS BEING STARTED OR NOT.)- 

■C THE VARIABLES WHICH MUST BE SAVED ARE 8 

C NE *NV* KDS *KDMAX *KSOUT * LDOU 8 *LFD , LSC * LSTC , I FL , I f LS * KQM (INTEGERS ) . } 

G. * ERNDyQDEC * ERRMX-* E2HAVE * E2HF AC* E2HMA X , RNDC.. (REAL) ' 

C HH *TOUT * TL (DOUBLE PRECISION) 

C IF THE GST 0 P FEATURE IS USED* THE EVALUATION OF G REQUIRES AN 

C INTEGRATION* AND THIS INTEGRATION MAY RESULT IN ANOTHER GSTOP, 

C THEN SOME ADDITIONAL VARI ABLES MUST BE SAVED. 

. C IN MANY APPLICATIONS NE NE Q ) * NV 1 = SUM OF ORDERS OF THE 

C DIFFERENTIAL EQUATIONS) *KD$(=KD) ,AND KDMAX ( -MAX I MUM ORDER OF ANY 

*C DIFFERENTIAL EQUATION) WILL BE THE SAME FOR EVERY INTEGRATION* 

-C AND-- HENCE NEED NOT BE SAVED. ' 

C - 

* C INSTRUCTIONS FOR MAKING CERTAIN CHANGES IN THIS SUBROUTINE ARE 

C GIVEN THROUGHOUT THE LISTING. TO FIND THESE INSTRUCTIONS* SEE 

C BELOW. • 

C 

C TO ELIMINATE THE GSTOP CAPABILITY, SEE JUST BELOW CARDS SEQUENCED.. 

C 541,703*745,791,1201*1234,1239,1253, AND 1330. 

C THIS MAKES THE SUBROUTINE SHORTER AND REDUCES OVERHEAD A LITTLE. 

-C - • 

C TO REMOVE THE INTERPOLATION CAPABILITY* SEE JUST BELOW CARDS 

C SEQUENCED 1070 AND 1219. . . • 

C- - .. .THE GSTOP FEATURE- MUST. ALSO. BE ELI MI NA TED.. SI NCE IT REQUIRES THE ..... 

C INTERPOLATION CAPABILITY* IF. OUTPUT POINTS ARE NOT HIT EXACTLY 

C (THEY ARE HIT EXACTLY IF HM AXA . LE . AB S ( DELT ) , AND INITIAL H= 

. C DELT*{ 2**(-K> ) , K=0 ,1,2...), THEN I F L AG=3 - ON THE FIRST STEP THAT 

C (T-TOUT )*H.GT.O (SEE THE USAGE OF DELT). I F LAG IS SET EQUAL TO -4 

C * ON THE LAST STEP THAT ( T-TFI NAL ) *H.LE.O . 

fC 

C THE OUTPUT OPTION GIVES OUTPUT OF VARIABLES USED IN THE 

C INTEGRATION ON EVERY STEP THAT NEQ.LE.O. (WHICH OF COURSE MUST 

•C BE SET AFTER THE INITIAL CALL TO THE INTEGRATOR) TO ELIMINATE 

C THIS OPTION, SEE JUST BELOW CARDS S EQU ENCED ;834 AND 1057. 

C 

C THE CHECK OPTION WHEN ADDED TO THE OUTPUT OPTION OUTPUTS EVERY 

C VARIABLE IN THE CALLING SEQUENCE JUST AFTER ENTERING AND JUST 

C BEFORE LEAVING THE INTEGRATOR WHEN NEQ=0 «. THIS OUTPUT IS 

C SOMETIMES USEFUL IN DEBUGGING A PROGRAM. TO INCLUDE THIS OPTION 

C SEE JUST BELOW CARDS SEQUENCED 613 AND 802. 
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